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1.  INTRODUCTION 

1.1  The  Skyshine  Problem 

Gamma-ray  doses  outside  of  areas  containing  nuclear  materials  can  generally 
be  broken  into  two  components.  The  first  component  is  the  direct  dose 
contribution  arising  from  gamma  photons  that  travel  directly  from  the  source  to 
the  detector  location.  The  direct  component  can  usually  be  evaluated  easily  using 
ray  analysis  techniques  [Ch84].  The  second  dose  component,  which  is  generally 
more  difficult  to  calculate,  is  due  to  indirectly  scattered  radiation  and  includes 
skyshine  radiation,  radiation  streaming  through  ducts,  and  radiation  reflected  from 
surfaces  (albedo  radiation)  [Ch84].  Skyshine  dose  refers  to  the  dose  caused  by  the 
reflection  of  photons  in  the  air  back  to  a  ground  target.  In  this  study,  approximate 
methods  for  calculating  the  indirect  skyshine  gamma  dose  are  considered. 

Outside  of  nuclear  facilities,  the  skyshine  dose  can  become  an  important 
component  of  the  total  offsite  dose  rate  [Pe84,  An87]  and  has  become  an  important 
concern  in  the  radiological  assessment  of  these  facilities.  Skyshine  dose 
calculations  are  required  for  accident  analysis  calculations  at  PWR  and  BWR 
nuclear  power  plants  [Pe84,  An87].  In  BWR  power  plants  the  movement  of 
nitrogen-16  from  the  reactor  to  the  turbine  room  also  requires  an  estimate  of  the 
skyshine  dose  to  be  made  during  normal  operations  [An87].  The  storage  of  nuclear 
waste  both  above  ground,  in  buildings,  and  below  ground  will  likewise  require  an 
analysis  of  the  skyshine  dose  during  the  design  of  the  disposal  facility. 

Evaluation  of  the  skyshine  dose  is,  in  general,  more  complicated  than  the 
ray-analysis  techniques  used  for  analyzing  the  direct  dose  [Fa87].    The  techniques 


used  to  evaluate  the  skyshine  dose  can  range  from  engineering  approximations 
[Pe84]  to  complete  numerical  solutions  of  the  multigroup  transport  equation. 

1.2  Previous  Skyshine  Investigations 

Using  a  Monte-Carlo  method,  Lynch  et  al.  [Ly58]  calculated  (at  various 
source-to-detector  distances)  the  dose  rate  caused  by  multiply  scattered  gamma 
photons  emitted  by  a  monoenergetic  monodirectional  beam  of  gamma  photons  into 
an  infinite  air  medium.  The  dose  rates  caused  by  photons  of  a  given  direction  and 
energy  were  normalized  to  unit  source  strength  (1  photon  per  unit  time)  to  yield  a 
response  function  for  a  source  photon  of  a  given  energy  and  direction  of  travel. 
Lynch' s  Monte—Carlo  calculations  were  carried  out  for  specific  photon  energies, 
specific  source-detector  distances,  and  specific  directions  of  travel  relative  to  the 
source-detector  axis  [Ly58].  The  response  functions  generated  from  the  Monte 
Carlo  results  are  valid  for  source-to-detector  distances  from  5  m  to  100  m,  for 
energies  from  0.6  MeV  to  12  MeV,  and  for  beam  angles  from  1  to  180  degrees 
relative  to  the  source-detector  axis  [Ly58]. 

Trubey  [Tr61]  proposed  an  approximate  method  for  calculating  the  skyshine 
flux  by  considering  only  a  single  scatter  for  the  gamma  photons.  This 
single-scatter  approximation  ignored  the  contribution  of  multiply  scattered  gamma 
photons  as  well  as  the  attenuation  and  buildup  of  photons  in  the  air.  Moreover,  it 
was  not  considered  suitable  for  a  shielded  source.  However,  the  single  scatter 
approximation  for  a  bare  source  agreed  well  with  the  Monte-Carlo  results 
calculated  by  Lynch  [Tr61]. 

Kitazume  [Ki68]  later  extended  the  single  scattering  approximation  by 
incorporating    exponential    attenuation    and    a    Taylor-type    buildup    factor. 


Exponential  attenuation  was  included  along  both  the  uncollided  photon  path  and 
the  scattered  photon  path.  The  Taylor  buildup  factor  was  applied  only  along  the 
scattered  photon  path.  The  inclusion  of  exponential  attenuation  and  photon 
buildup  in  the  single  scattering  formulation  allowed  this  method  to  be  used  beyond 
the  100m  source—  to-detector  lengths  considered  by  Trubey  [Ki68].  Kitazume's 
results  were  in  good  agreement  with  Lynch's  Monte  Carlo  results  for  beam  angles 
less  than  60  degrees  at  all  energies  and  source—to-detector  distances  considered. 
At  angles  greater  than  60  degrees,  Kitazume's  method  was  in  better  agreement 
with  Lynch's  results  than  were  Trubey's  results,  being  at  most  20%  lower  than 
results  reported  by  Lynch  et  al.  [Ki68]. 

With  these  successful  point-kernel  approaches  for  estimating  the  photon 
skyshine  dose,  general  purpose  codes  were  developed  for  use  in  complex  physical 
geometries.  One  such  widely  used  point-kernel  code  is  G3  [Ma73]  which  uses 
surfaces  defined  by  quadratic  equations  to  define  the  geometry  of  the  problem.  G3 
uses  exponential  attenuation  of  the  direct  (uncollided)  beam  and  has  the  option  of 
using  buildup  of  scattered  photons  in  the  scattered  leg  to  account  for  multiple 
scattering  [Ma73j.  However  at  present,  G3  cannot  deal  with  buildup  in  shielding 
structures  along  the  unscattered  photon's  path  of  travel. 

A  specialized  extension  of  the  single-scattering  approximation  accounting  for 
overhead  shielding  above  a  point  isotropic  source  was  proposed  by  Roseberry  and 
Shultis  [Ro80,  Ro82].  Roseberry  and  Shultis  proposed  a  point-kernel  model  with 
exponential  direct  beam  attenuation,  buildup  of  scattered  photons  in  the  scattered 
air  leg,  and  an  infinite-medium  buildup  factor  for  the  overhead  shield.  The 
infinite-medium  buildup  factor  was  introduced  to  approximate  the  scattered 
photons  produced  in  the  shield.     Their  model  gave  reasonably  good  agreement 


when  compared  to  experimental  data  obtained  from  a  shielded  skyshine  benchmark 
experiment  [C178,  Sh78,  Na81]. 

A  more  accurate  way  of  calculating  the  skyshine  dose  for  an  overhead 
shielded  gamma  source  would  be  to  use  a  calculational  technique  based  on  an  exact 
transport  description  of  the  particular  skyshine  problem.  The  solution  of  the 
transport  equation  for  skyshine  geometries  normally  requires  a  multidimensional 
geometric  representation  and  a  multigroup  energy  formulation.  One  way  for 
solving  such  a  multidimensional,  multigroup  transport  equation  is  to  use  Monte 
Carlo  techniques.  General  purpose  Monte  Carlo  codes  that  are  available  to  solve 
skyshine  problems  include  COHORT  [So75]  and  OGRE  [Pe65]. 

Discrete  ordinates  codes  such  as  DOT  [My73]  and  ANISN  [En67]  have  also 
been  used  to  solve  some  skyshine  problems.  The  discrete  ordinates  procedure  can 
give  very  accurate  solutions  for  a  geometry  (usually  simple)  that  it  is  capable  of 
modeling.  Unfortunately,  both  multidimensional  discrete-ordinates  solutions  and 
multidimensional  Monte-Carlo  results  require  large  computational  effort,  thereby 
limiting  the  usefulness  of  these  codes  for  preliminary  or  routine  design  and  safety 
studies. 

To  reduce  the  computational  effort,  and  yet  maintain  acceptable  accuracy  for 
the  estimation  of  skyshine  doses  from  shielded  sources,  codes  based  on  semi- 
empirical  skyshine  methods  have  been  developed.  Such  specifically  designed 
skyshine  codes  use  response  functions  (obtained  by  fitting  an  empirical  formula  to 
Monte  Carlo  skyshine  results)  to  calculate  the  skyshine  dose.  Examples  of  such 
codes  include  SKYSHINE  [Pr76],  SKYSHINE-II  [La79],  and  SKYSHINE-III 
[La88].  These  three  codes  consider  radiation  sources  (photons  or  neutrons)  in  a 
rectangular  structure  with  four  walls,  a  roof,  and  a  floor.   The  codes  break  each  of 


the  containment  surfaces  up  into  a  series  of  different  sections,  each  with  its  own 
attenuation  properties.  With  a  Monte  Carlo  sampling  technique,  the  radiation 
energy  and  location  on  the  containment  surface  through  which  the  radiation  will 
stream  are  chosen.  After  a  correction  for  attenuation  as  the  beam  penetrates 
through  the  structure,  the  contribution  to  the  skyshine  dose  made  by  the 
transmitted  beam  is  then  calculated  with  the  use  of  the  beam  response  formulas. 

To  reduce  the  computational  effort  of  running  multidimensional  transport 
codes  for  shielded  skyshine  problems,  a  one-dimensional  transport  code  can  be  run 
to  determine  the  photon  distribution  leaving  the  shielding  around  a  source.  This 
source  distribution  can  then  be  used  with  any  skyshine  calculational  method 
suitable  for  treating  the  unshielded  skyshine  problem.  This  hybrid  approach  was 
used  by  Keck  and  Herchenroder  [Ke82]  who  combined  the  ANISN  and  the 
SKYSHINE  II  codes  to  calculate  skyshine  dose  for  the  K-State  Benchmark 
experiment  [C178].  Peng  [Pe84]  also  followed  this  two-stage  approach  by  using 
ANISN  and  COHORT  to  calculate  the  skyshine  dose  rate  outside  of  a  nuclear 
power  plant  during  a  LOCA  accident  analysis. 

To  reduce  the  cost  of  analyzing  shielded  skyshine  doses,  Faw  and  Shultis 
[Fa87,  Sh87]  recently  modified  the  original  SKYSHINE-II  method  [La79].  They 
developed  improved  beam  response-function  formulas  by  fitting  a  three-parameter 
formula  to  skyshine  results  obtained  with  a  point-kernel  technique  which 
accounted  for  gamma-ray  attenuation,  photon  pair  production,  and  buildup  in  the 
scattered  leg.  Unlike  the  SKYSHINE-II  method  (which  used  Monte  Carlo 
techniques  to  account  for  different  source  emission  directions),  the  skyshine  dose 
was  found  by  numerically  integrating  over  all  emission  directions.  This  revised 
skyshine  method   was  incorporated   in   the   microcomputer   code   MicroSkyshine 


[Gr87].  The  MicroSkyshine  method  also  added  an  interpolation  scheme  to  make 
the  beam  response  functions  continuously  variable  in  both  energy  and  angle.  The 
new  response  functions  also  eliminated  the  stochastic  variations  observed  in  the 
response  functions  used  in  the  original  SKYSHINE-II  method. 

MicroSkyshine  can  treat  point,  isotropic,  polyenergetic,  gamma  sources  with 
or  without  overhead  shielding  in  two  basic  geometries.  The  first  geometry 
MicroSkyshine  treats  is  the  case  of  a  gamma  photon  source  located  on  the  axis  of  a 
cylindrical  silo.  The  second  geometry  has  the  source  and  detector  located  on 
opposite  sides  of  a  vertical  wall  which  may  be  oriented  obliquely  to  the 
source-detector  axis. 

MicroSkyshine  calculates  the  skyshine  dose  by  integrating  the  fitted  beam 
response  functions  over  all  directions  allowed  by  the  geometry  of  the  problem.  The 
effect  of  any  overhead  shielding  is  accounted  for  by  exponentially  attenuating  the 
beam  through  the  shield  and  then  correcting  the  attenuation  in  the  shield  by 
multiplying  by  a  buildup  factor  to  obtain  an  estimate  of  all  the  gamma-rays 
(uncollided  or  scattered)  passing  through  the  shield.  The  MicroSkyshine  method 
was  found  to  be  in  good  agreement  with  other  methods  and  with  the  K-State 
benchmark  skyshine  experiment  [C178]  for  the  unshielded  cases. 

For  shielded  skyshine  problems,  there  is  a  sparsity  of  published  calculations 
and  skyshine  experimental  measurements.  The  K-State  benchmark  skyshine 
experiment  included  two  shielded  silo  configurations.  Although  the  MicroSkyshine 
code  gave  better  agreement  with  these  experimental  results  than  did  other 
calculational  methods  [Fa87j,  the  accuracy  and  robustness  of  the  MicroSkyshine 
method  for  treating  sources  of  different  energies  and  degrees  of  collimation  and 
shields  of  different  materials  and  thicknesses  was  largely  uncertain. 


1.3  Purpose  of  Study 

This  study  was  motivated  by  the  need  to  investigate  the  capabilities  and 
limitations  of  the  MicroSkyshine  method  for  treating  skyshine  sources  with  a  slab 
overhead  shield.  Inherently  accurate  methods  for  calculating  the  skyshine  dose  are 
methods  based  on  a  complete  multidimensional  description  of  the  radiation  field 
using  the  photon  transport  equation.  However,  the  large  computational  cost 
needed  to  solve  numerically  the  multidimensional  transport  equation  precludes 
using  this  approach  to  obtain  the  many  benchmark  calculations  needed  to  assess 
the  MicroSkyshine  method. 

Less  expensive  approaches  for  calculating  the  skyshine  dose  from  a  shielded 
source  include  (1)  the  buildup  factor  approach  for  the  shield  used  by  Roseberry 
[R08O,  Ro82]  and  Faw  and  Shultis  [Fa87,  Sh87],  and  (2)  the  composite  method 
used  by  Keck  and  Herchenroder  [Ke82]  and  Peng  [Pe84]. 

Although,  the  buildup  factor  approach  is  the  simplest  approach  to  the 
shielded  skyshine  problems,  the  magnitude  of  the  error  introduced  by  this  method's 
approximations  remains  unexplored  and  unverified.  Consequently,  in  this  study  a 
composite  method,  similar  to  Keck  and  Herchenroder' s  [Ke82],  was  developed. 
The  composite  method  uses  an  accurate  one-dimensional  transport  code  to 
compute  the  emergent  energy-angle  distribution  of  photons  on  the  outer  shield's 
surface.  Then  a  modified  MicroSkyshine  method  uses  this  emergent  distribution  of 
photons  as  a  bare  skyshine  source  to  calculate  the  skyshine  dose  at  the  detector 
location.  This  composite  method,  once  verified,  is  used  to  investigate  the  accuracy 
of  the  buildup-factor  approach  used  by  MicroSkyshine  to  estimate  the  skyshine 
dose  arising  from  a  shielded  point  source  in  a  silo. 


The  MicroSkyshine  method  for  calculating  the  skyshine  dose  is  reviewed  in 
Chapter  2.  The  discussion  of  the  MicroSkyshine  method  focuses  on  the  generation 
of  the  fitted  beam  response  functions,  the  energy  and  angular  interpolation  of  the 
response  functions,  and  the  process  by  which  the  response  functions  are  integrated 
to  obtain  the  skyshine  dose  at  the  detector.  In  addition,  an  extension  of  the 
original  MicroSkyshine  method  to  treat  sources  with  variable  angular  intensities  is 
presented  in  Chapter  2. 

In  Chapter  3,  the  concept  of  using  a  two-step  composite  method  to  calculate 
the  skyshine  dose  is  explored.  Specifically,  this  composite  method  combines  a 
one-dimensional,  discrete-ordinates,  transport  equation  for  the  source  shield  with 
an  extension  of  the  MicroSkyshine  method  for  calculating  the  skyshine  dose  from 
unshielded  sources.  The  validity  of  using  this  two-step  scheme  for  calculating  the 
skyshine  dose  is  explored  and  shown  to  be  very  reasonable  for  the  class  of  skyshine 
problems  considered  in  this  study. 

Also,  in  Chapter  3,  a  procedure  to  allow  the  use  of  a  one-dimensional 
transport  equation  for  the  inherently  two-dimensional  transport  problem  of  a  point 
source  shielded  by  a  cylindrical  slab  is  presented.  This  transformation  technique 
allows  the  use  of  a  one-dimensional,  discrete-ordinates,  transport  calculation, 
rather  than  a  two-dimensional  calculation,  to  estimate  the  source  distribution 
emerging  from  the  slab  shield,  i.e.  the  distribution  needed  by  the  modified 
MicroSkyshine  method.  The  basic  outline  of  the  discrete  ordinates  method  for 
solving  a  one-dimensional  transport  problem  is  reviewed.  The  needed  boundary 
source  term  for  the  one-dimensional  discrete  ordinates  transport  solution  is  then 
derived  for  the  case  of  slab  shielding  over  the  skyshine  source,  and  the  process  of 
linking  the  output  from  the  discrete  ordinates  code  to  the  modified  MicroSkyshine 


method  is  explained.  Finally  in  Chapter  3,  comparisons  are  given  between  the 
composite  method  to  results  from  the  K-State  Benchmark  Experiment. 
Comparisons  to  other  skyshine  calculation  methods  are  also  presented. 

In  Chapter  4,  the  composite  method  is  used  to  investigate  the  accuracy  of  the 
MicroSkyshine  method  (exponential  attenuation  with  buildup  in  the  shield)  for 
estimating  the  skyshine  dose  from  a  shielded  point  source.  In  this  investigation 
gamma  ray  energies  of  6.129  MeV,  1.25  MeV,  and  0.5  MeV,  and  4  different  silo 
geometries  with  various  shield  thicknesses  are  considered. 

Finally,  Chapter  5  presents  conclusions  reached  during  this  investigation  as 
well  as  suggestions  for  further  study. 


2.   REVIEW  AND  MODIFICATION  OF  THE  MICROSKYSHINE  METHOD 

The  beam  response  functions  used  by  Faw  and  Shultis  [Fa87,  Sh87]  to 
calculate  the  skyshine  dose  in  the  MicroSkyshine  method  where  generated  by 
fitting  an  empirical,  three-parameter,  response-function  formula  to  results 
obtained  by  applying  the  point  kernel  method  to  a  monoenergetic  monodirectional 
beam  of  photons.  Figure  2.1  shows  the  geometry  used  in  estimating  the  response 
at  a  point  isotropic  detector  caused  by  the  monoenergetic  monodirectional  beam  of 
photons  in  an  infinite,  homogeneous,  air  environment.  The  response  produced  by 
the  detector  will  depend  upon  the  photon's  initial  energy,  the  photon's  initial 
direction  of  travel  with  respect  to  source-detector  axis,  the  material  in  which  the 
photon  travels,  and  the  source-to-detector  distance. 

In  this  chapter,  the  MicroSkyshine  methodology  for  calculating  the  skyshine 
dose  is  studied.  In  particular,  the  review  of  the  MicroSkyshine  method  examines 
the  generation  of  the  beam  response  functions,  the  fitting  of  an  empirical  formula 
to  these  response  functions,  and  the  use  of  the  response  functions  in  calculating  the 
skyshine  dose.  The  MicroSkyshine  method  is  then  extended  to  treat  anisotropic 
point  sources  with  variable  source  energies. 

2.1  Calculation  of  Detector  Response  to  a  Gamma  Photon  Beam 

The  probability  a  photon  (see  Fig.  2.1)  of  energy  E  travels  a  distance  y  in  air 
without  interaction  and  then,  while  traveling  a  further  distance  dy,  scatters 
through  an  angle  9S  into  a  unit  solid  angle  is  [Ch84] 

Z  N  e<Tc(E,0s)  e~^y  dy,  (2.1) 
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Figure  2.1.    Geometrical  representation  used  in  calculating  the  response  functions 
for  MicroSkyshine  [Sh87]. 
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where  ZN  is  the  electron  density  of  the  air,  //  is  the  total  interaction  coefficient  at 
energy  E,  and  e0c(E,0s)  is  the  Klein-Nishina  differential  scattering  cross  section. 
This  differential  scattering  cross  section  can  be  approximated  by  the  free-electron 
Klein-Nishina  differential  scattering  cross  section,  and  is  given  by  [Ch84] 

e<7c(E,0s)  =  \  re2  p  [1  +  p2  -  p  (1-COS0S)],  (2-2) 


where 


n  _  mec2  /.-)  o\ 

p-mPc2  +  E  -  Ecos0s  '  (2-,3) 


mec2  =  0.511  MeV,  and  r  is  the  classical  electron  radius.  Source  photons  with 
energies  greater  than  1.02  MeV  can  produce  annihilation  photons.  The  probability 
that  a  source  photon  while  traveling  a  distance  dy  after  traveling  a  distance  y 
produces  an  annihilation  photon  traveling  towards  the  detector  is  [Sh87] 


^  N  (7pp(E)  e  W  dy  ,  (2.4) 


where  N  is  the  atomic  density  of  the  material  being  travelled  through  and  crPP(E) 
is  the  microscopic  pair  production  cross  section  at  energy  E. 

The  response  at  the  detector,  without  buildup,  due  to  a  beam  of  photons 
with  energy  E  and  with  direction  of  travel  0  can  be  found  by  integrating  Eqs.  (2.2) 
and  (2.4)  along  the  path  of  the  photon  beam  to  obtain 
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.00    -fjcy  , 

5&(E,0,x)=Nf    VPe^WRtEV 
J0      r 


+  ^<7PP(E)e^rR(Ea)]dy.  (2.5) 


The  scattered  photon  energy  is  denoted  by  E',  the  photon  annihilation  energy  by 
Ea,  the  scattering  point  to  detector  distance  by  r  and,  the  detector  "response 
function"  at  a  scattered  energy  Ed  by  R(E<i).  The  scattered  energy  E'  is  found 
using  the  Compton  formula  for  free  electron  scattering  [Ch84],  namely, 

p 
E'  =  1  +  (E/mec2)(l  -  cos^s)  •  (2-6) 

To  account  for  subsequent  buildup  of  scattered  photons  along  the  second  (or 
scattered)  leg,  the  scattering  source  term  and  the  pair  production  term  in  Eq.  (2.5) 
was  multiplied  by  an  appropriate  buildup  factor  B(E,r).  Thus,  the  total  detector 
response  is  estimated  by 


oo    —fjcy  i 

a(E,0,d)  =  N/    ^  [Ze<7c(E,0s)R(E')B(E',r)e  ^r 

Jo     r 


+  B(Ea^ri(7pp(E)R(Ea)e-//ar]dy>  (2>7) 


2tt 


A  change  in  the  variables  of  integration  using  y  =  /xy,  r'  =  /i'r,  r"  =  //ar,  and  -z 
lii  results  in  a  detector  response  of,  [Sh87] 
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oo      -y 
tf(E,0,d)  =  N  /    =V[Ze(7c(E^s)B(E',r') 


+  B(E^lR(Ea)£rpp(E)e-r]dy  (2>g) 


To  evaluate  Eq.  (2.8),  the  buildup  factor  and  interaction  coefficients  must  first  be 
selected.  The  total  mass  interaction  coefficients  were  taken  from  Hubbell  [Hu82] 
and  the  annihilation  microscopic  cross  section  data  were  taken  from  Storm  and 
Israel  [St67].  The  buildup  factor  B(E,/zr)  was  assumed  equal  to  the  infinite 
medium  exposure  buildup  factor  as  approximated  by  the  geometric  progression 
model  [Ha83,  Ha86]  for  a  point  isotropic  source.  The  geometric  buildup-factor 
model,  used  in  evaluating  the  beam  response  functions,  is  given  by 


B(E,X)  = 


1  4-  (b-l)(KX-l)  ,  for  k  }   1 
.1  +  (b-1)  X  ,  fork=l 


(2.9) 


where  X  =  fix  is  the  source-to-detector  distance  in  mean-free-path  (mfp)  lengths, 
and  K  is  evaluated  from 

K(X)  =cXa  +  d  ^^/^tShjl^f  ("2)1  •  (2-10) 

Values  for  the  coefficients  a,  b,  c,  d,  and  Xk  were  taken  from  a  recent  revision  of 
QAD  [Rs86]  and  depend  on  the  photon  energy  and  the  shielding  material  [Sh87]. 

Evaluation  of  the  integral  in  Eq.  (2.8),  once  the  integrand  was  known,  was 
performed  using  Gaussian  quadrature  [Sh87].  A  16-point  Gaussian  quadrature 
routine  integrated  the  detector  response  for  each  mfp  of  interest  along  the  beam 
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path  of  the  source  photon.  The  integration  procedure  started  at  the  source  and 
continued  along  the  photon  path-of-travel  until  the  change  in  the  value  of  the 
integral  was  less  than  a  prescribed  value  [Sh87].  In  this  manner  the  detector 
response  function  was  evaluated  at  a  given  set  of  energies  and  for  each  energy  at  a 
given  set  of  angular  directions.  Table  2.1  shows  the  energy  set  and  Table  2.2 
shows  the  angular  set  used  in  evaluating  the  beam  response  functions  [Sh87]. 

Several  important  approximations  were  made  in  deriving  Eq.  (2.7)  for  the 
beam  response  functions  [Sh87].  First,  the  Klein-Nishina  differential  scattering 
cross  section  ignores  all  electron  binding  effects  on  the  cross  sections  [Ch87]. 
Errors  caused  by  the  Klein-Nishina  approximation  over  the  energy  range  of 
interest  (0.1  <  E  <  10  MeV)  are  expected  to  be  very  small  since  electron  binding 
effects  are  relatively  small  at  these  relatively  high  energies.  Second,  the  buildup 
factors  used  were  based  on  an  isotropic  point  source,  a  situation  not  rigorously 
realized  for  this  analysis,  since  photons  which  scatter  in  dy  are  preferentially 
scattered  in  a  forward  direction.  Thus,  use  of  the  buildup  factor  for  a  point 
isotropic  source  will  tend  to  overestimate  the  dose  at  the  detector. 

However,  the  errors  introduced  by  the  model  assumptions  discussed  above 
appear  to  be  quite  small  because  of  the  excellent  agreement  between  benchmark 
experimental  data  and  the  skyshine  calculations  using  the  above  beam  response 
functions  [Sh87]. 
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Table  2.1.  Energy  group  structure  used  in  deriving 
the  beam  response  functions  used  in 
MicroSkyshine  [Sh87] . 


Energy  Group 

Group  Energy 
(MeV) 

1 

1 

9.5 

2 

8.5 

3 

7.5 

4 

6.5 

5 

5.5 

6 

4.5 

7 

3.5 

8 

2.5 

9 

1.5 

10 

0.75 

11 

0.325 

12 

0.055 

Table  2.2    The  discrete  beam  directions  used  by  Faw  and 
Shultis  in  deriving  the  new  beam  response 
functions  used  in  MicroSkyshine  [Sh87] . 


Angular 

Angle 

Angular 

Angle 

Group 

<t>) 

Group 

0j 

1 

0.5 

11 

45.0 

2 

1.5 

12 

55.0 

3 

2.5 

13 

65.0 

4 

4.0 

14 

75.0 

5 

6.0 

15 

85.0 

6 

8.5 

16 

95.0 

7 

12.5 

17 

110.0 

8 

17.5 

18 

130.0 

9 

25.0 

19 

150.0 

10 

35.0 

20 

170.0 
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A  third  important  approximation  in  the  MicroSkyshine  method  is  the  use  of 
an  infinite  air  medium  for  evaluating  the  beam  responses.  A  more  correct 
approach  would  include  the  effect  of  the  ground  in  the  beam  response  calculations. 
However,  the  complications  involved  in  adding  an  air-ground  interface  are 
considerable  and  include  changes  in  the  detector  response  for  different  detector 
heights  above  the  ground  and  different  soil  compositions.  However,  since  the 
average  Z  number  for  most  soils  is  reasonably  close  to  that  of  air,  and  since  photon 
reflection  from  the  ground  is  typically  small  compared  to  the  contribution  of 
photons  approaching  the  detector  from  above,  the  neglect  of  the  air-ground 
interface  is  felt  to  be  justified  [Sh87] . 

2.2  Approximation  of  Beam  Response  Functions 

To  simplify  the  use  of  beam  response  functions,  Faw  and  Shultis  [Fa87,  ShS7] 
fit  a  semi-empirical  formula  to  their  calculated  values  of  the  detector  response. 
The  detector  responses  calculated  by  the  point  kernel  were  approximated  using 
[Sh87] 

#(E,^,x)  =  E  <?(E,0,x)  ,  (2.11) 

where  the  fit  formula  was 

^(E,0,x)  =  Mp/Pof  [x  (p/p0)]h  e(a^X  p/po).  (2.12) 

The  fit  parameters  a,  b,  and  c  depend  only  upon  source  energy  E  and  beam 
direction  0.  The  reference  air  density  is  given  by  p0  and  the  air  density  to  be  used 
in  a  dose  calculation  by  p.    The  reference  air  density  was  0.001225  g/cm3  [Sh87]. 
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The  detector  response  #(E,0,x)  as  calculated  has  units  of  rad  per  photon  when 
X-  1.308X10"11  rad  m2/MeV  [Sh87]. 

Values  of  the  fit  parameters  a,  b,  and  c  were  found  for  a  fixed  Ej  and  0j  by 
minimizing  the  squared  sum  of  the  difference  between  the  calculated  response 
functions  and  the  prediction  from  Eq.  (2.11)  for  all  calculated  source-to-detector 
distances  [Sh87].  Since,  the  response  at  the  detector  was  found  to  vary  over  many 
orders  of  magnitude  as  the  source-to-detector  distance  x  changed,  the  least  squares 
fits  were  actually  performed  by  fitting  the  natural  log  of  Eq.  (2.11)  to  the  natural 
log  of  the  calculated  responses  [Sh87].  Using  the  logarithmic  fitting  procedure  for  p 
=  p0  results  in  the  following  least  squares  fit  function. 

M 

Sij(a,b,c)  =      I   [G  +  b  ln(Xm)  +  a  -  cxm  -  /rc{Rm(Ei,0j,xm)}p,         (2.13) 
m=l 

where  G  =ln(E{)  and  -5£(Ei,0j,xn)  is  the  beam  response  at  a  detector  distance  of 
xm.  The  parameters  a,  b,  and  c  that  minimized  S  were  then  found  using  the 
simplex  method  [Sh87]. 

The  point-kernel  beam  responses  were  calculated  out  to  source-to-detector 
distances  of  2500m.  Then,  the  fitted  response  parameters  were  calculated  for  all 
energy  and  angular  values  given  in  Tables  2.1  and  2.2.  Comparison  of  the  fitted 
(approximate)  response  functions  with  the  calculated  values  indicated  that  the 
poorest  agreement  always  occurred  at  small  source-to-detector  distances.  The 
error  associated  with  the  fitted  response  functions  was,  for  almost  all  the  cases,  less 
than  an  absolute  deviation  of  16%.  The  detailed  results  calculated  by  Faw  and 
Shultis  are  in  reference  works  by  those  authors  [Sh87]. 


18 


2.2.1   Interpolation  of  the  Fitted  Response  Functions 

An  interpolation  scheme  was  used  in  MicroSkyshine  to  evaluate  the  beam 
response  functions  at  energies  and  angular  directions  different  from  those  used  to 
approximate  the  response  functions  [Sh87].  The  interpolation  scheme  results  in 
beam  response  functions  that  are  continuously  variable  in  energy  and  direction. 

The  beam  response  function  for  a  photon  with  energy  E  and  direction  of 
travel  (f>  is  found  using  a  linear  interpolation  scheme.  The  beam  response  function 
is  first  interpolated  in  energy  at  all  discrete  angular  directions  0j  using  [Sh87] 


^(E,^,x)  =  ^i+1J 


Ej  -  E 
Ei-  E 


i+1 


Mj 


E"Ei  +  l 
Ei  -  E 


i+1 


(2.14) 


where 


and 


Jf.  =  ^(Ei,0j,x)  , 


(2.15) 


Jf+1J  ~=  '(Ej+pty*)  •  (2-16) 

For  energies  between  9.5  and  10.0  MeV,  the  following  extrapolation  scheme  is  used 
[Sh87] 


^(E,0j,x)  =  X  •  (E  -  8.5)  +  3L  .  (9.5  -  E) 


(2.17) 


Then,  the  response  function  for  photons  of  energy  E  with  direction   0  (where 
0j  <  (j>  <  <j>.    A  is  estimated,  using  linear  interpolation,  as  [Sh87] 

^(E,0.   15x)  -  ^(E,0j,x) 
^(E,0,x)  =  ^(E,0j,x)  + )+\        _  6 [<t>-<P-]+1\-      (2-18) 
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For  directions  in  the  two  end  intervals,  linear  extrapolation  is  used,  namely,  for 
1700  <  0  <  1800, 

*(E,*,x)  =  ,?(E,02O,x)  +  *l^fi)9  I  fflfr»i«)  (4>-M  ■       (2.19) 

and  0°  <  0  <  .  50 

^(E,0,x)  =  <T(E,fr,x)+^A^  I  |SAi*)  (0-02)  .  (2.20) 

2.3  Skyshine  Dose  Calculations  Using  Beam  Response  Functions 

The  skyshine  dose  at  a  source-to-detector  distance  of  d  arising  from  a 
collimated,  bare,  isotropic,  point  source  can  be  found  by  integrating  the  response 
functions  over  all  source  energies  and  over  all  possible  photon  beam  directions  Qs. 
The  skyshine  dose  can  thus,  be  calculated  from 


.00 


R(d)  =  f    dE'  f     dQ  S(E',n)  5&(E',0,d)  ,  (2.21) 

■'o       J  ns 


where  S(E',fl)  is  the  energy  and  directional  distribution  of  the  point  skyshine 
source,  Qs  is  the  source's  solid  angle  of  collimation,  and  ^(E',0,d)  is  the  beam 
response  function  which  may  be  approximated  by  Eq.  (2.11).  Often,  the  source 
S(E,fi)  is  a  monoenergetic  isotropic  point  source  with  energy  E.  For  this  case 
S(E\fi)  becomes 


S(E',fi)  =  |f  <5(E'-E)  ,  (2.22) 
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and  Eq.  (2.21)  reduces  to 


R(d)  =  |f  /     m,M  <M  •  (2.23) 

Sis 


2.3.1  Skyshine  Calculations  for  a  Source  in  a  Silo 

The  geometry  of  the  silo  problem  is  illustrated  in  Fig.  2.2.     For  the  silo 
problem,  the  various  geometric  variables  are  defined  as 

ys  =  the  source  height  below  the  silo  top, 

yd  =  the  detector  height  below  the  silo  top, 

x  =  the  horizontal  source-to-detector  distance, 

r  =  the  silo  radius, 

^o  =  /  1  +  r2/y^  =  cos  60  =  the  source  collimation  angle, 


d     =  J  x'2  +  (yd-ys)2  =  the  source-to-detector  distance, 


and 


a    =  cos_1(x/d)  =  the  angle  between  the  horizontal  and  the 

source-detector  axis. 
The  solid  angle  of  collimation  for  the  source  is  given  by  2lj0.  Evaluation  of 
Eq.  (2.23)  for  the  Skyshine  dose  for  the  silo  problem  first  requires  that  the  angle  0 
be  expressed  in  terms  of  the  integration  variables.  From  Figure  2.3  it  is  seen  that 
cos0  is  the  dot  product  of  a  unit  vector  along  the  photon's  initial  path  and  the  unit 
vector  along  the  source—to-detector  axis.  The  angle  0  in  terms  of  6,  ip,  and  a  is 
thus  found  to  be 

0  =  cos'^sintf  cos^  coxa  -  cos#  sina)  .  (2.24) 
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Figure  2.2.    Illustration  of  the  variables  representing  the  physical  distances  used  in 
the  MicroSkyshine  method  for  a  point  source  in  a  silo  [Sh87]. 
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photon 
ray 


yd  _  ys 


detector 


Figure  2.3.    Angular  coordinate  system  used  to  transform  the  integral  equation  in 
MicroSky shine  [Sh87]. 
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2.3.2  Unshielded  Silo  Calculation 

The  skyshine  dose  for  an  open  silo  with  a  point  isotropic  source  is  found  from 
Eq.  (2.23).   For  the  silo  problem  the  differential  solid  angle  is  given  by 

dfl  =  sm6  d9  dip  .  (2.25) 

Substitution  of  Eq.  (2.25)  into  Eq.  (2.23)  results  in  the  detector  response  being 
given  by 


c       ,•  "max         /•^7I" 

R(d)  =  ±r  f        (±0  J       sin0  #(E,0,d)  dip  •  (2.26) 


In  the  silo  case,  the  azimuthal  contribution  is  symmetric  about  the  azimuthal 
reference  axis  (i.e.,  the  source-to-detector  axis).  This  azimuthal  angular 
symmetry  reduces  Eq.  (2.26)  to 


c       -  7T  -  #max 

R(d)  =  fc  J     dip  f         d$sm6  #(E,0,d)  ,  (2.27) 


or  if  one  lets  w  =  cos#  then  Eq.  (2.27)  reduces  to 


R(d)  =  §*  /    dip  f    du  #(E,0,d)  .  (2.28) 

0  LU0 
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Equation  (2.28)  may  be  evaluated  numerically  using  a  double  Gaussian  integration 
scheme  [Sh87]. 

2.3.3  Shielded  Silo  Dose  Calculation 

Faw  and  Shultis  [Fa87,  Sh87]  also  accounted  for  the  addition  of  a  horizontal 
slab  shield  on  top  of  the  silo  by  using  exponential  attenuation  and  a  buildup  factor 
along  the  path  through  the  shielding  slab.  The  path  length  through  the  shield  is 
illustrated  in  Fig.  2.4.  The  mean-free-path  distance  along  the  path  through  the 
shield  is 

where  t  is  the  shield  thickness,  ps  is  the  shield's  density,  (fi/p)  is  the  total  mass 
interaction  coefficient,  and  u  is  the  cosine  of  the  polar  angle  6.  Inclusion  of 
buildup  and  exponential  attenuation  in  the  shield  results  in  a  detector  response  of 


c       -7T  -  #max  _\ 

R(d)  =2£  f    df  f        do;B(e,A)  ^(E,0,d)  e  A  .  (2.30) 

Z7r*/0         ^0 


The  buildup  factor  B(E,A)  was  approximated  by  a  Berger  approximation  of  the 
form 

B(E,A)  =  1  +  aA  ebA  .  (2.31) 

Values  for  the  mass  interaction  coefficients  (ft/p)  for  air,  water,  concrete,  iron, 
lead,  zirconium,  and  uranium  dioxide  were  again  taken  from  data  by  Storm  and 
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Figure  2.4.    Effective  slab  thickness  for  a  photon  penetrating  a  slab  of  thickness  t 
at  an  angle  0. 
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Israel  [St67].  The  interpolation  in  energy  of  the  mass  interaction  coefficients  was 
performed  using  a  logarithmic  fit.  The  Berger  buildup  factor  B(E,A)  [Ch84]  is 
usable  for  energies  between  0.015  MeV  and  15  MeV  and  for  shield  thicknesses  up  to 
40  mfp.  The  buildup-factor  parameters  a  and  b  were  interpolated  in  energy  using 
a  linear  energy  fit. 

2.4  Extension  to  Polyenergetic  and  Anisotropic  Point  Sources 

The  original  MicroSkyshine  method  [Fa87,  Sh87]  was  limited  to  point, 
monoenergetic,  isotropic  sources.  To  use  the  MicroSkyshine  response  functions  for 
cases  involving  anisotropic  point  sources  with  variable  source  energies,  it  is 
necessary  to  modify  the  original  method.  This  section  presents  a  modification  for 
calculating  the  skyshine  dose  caused  by  an  anisotropic  and/or  variable  energy 
source. 

The  skyshine  dose  for  anisotropic  point  sources  with  a  distribution  of  source 
energies  is  formally  given  by  Eq.  (2.21).  When  an  unshielded,  anisotropic, 
polyenergetic,  point  source  is  placed  on  the  axis  of  the  cylinder,  the  detector 
response  thus  becomes 


-Emax  -27T  t/max 

R(d)=  f         dE  J       d^>  /  d0sin0S(E,0,^)  #(E,0,d)  ,  (2.32) 

^0  ^0  ^0 


where  0  is  given  by  Eq.  (2.24),  d  is  the  source-to-detector  distance,  and  S(E,0,^>) 
is  the  energy-direction  distribution  of  source  photons  emitted  with  energy  E  and  in 
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direction  defined  by  9  and  ip.   With  the  definition  u  =  cosO,  Eq.  (2.32)  reduces  to 


-Emax  /27T  1 

R(d)  =  J  dE  J       dip  J    du  S(E,w,tf>)  #(E,0,d)  .  (2.33) 

^0  ^0  Ju0 


Of  particular  interest  in  this  study  is  the  case  when  the  source's  angular 
distribution  is  independent  of  0  and  when  the  source's  continuous  energy 
distribution  is  approximated  by  a  multigroup  approximation.  The  multigroup 
approximation  of  the  source's  energy  distribution  can  be  incorporated  in  Eq.  (2.33) 
by  replacing  the  energy  integral  with  a  summation  over  the  midpoint  energy  of 
each  energy  group,  i.e., 

G         2j  1 

R(d)  =    I  dip  I     dw  Sg(u;)  #(Eg,0,d)  .  (2.34) 

g=l   J0  Ju0 

The  sources  angular  independence  of  ip  allows  Eq.  (2.34)  to  be  reduced  to 

G       *         1 

R(d)  =  2    J  <W        dwSg(u>)  #(Eg,<£,d)  ,  (2.35) 

g=l    0        Ju0 

when  dealing  with  an  azimuthally  symmetric  distribution  of  source  photons. 

The  numerical  evaluation  of  the  integrals  in  Eq.  (2.35)  can  be  done  using  a 
numerical  quadrature  method  such  as  Gaussian  quadrature.  Equation  (2.35)  can 
thus  be  approximated  by 
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G      N  N 

R(d)  =  2    I      I    wi    I    Wj  SgC^j)  R(Eg,0,d)  ,  (2.36) 

g=l  i=l       j=l 

where  wj  and  wj  are  the  Gaussian  quadrature  weights  and  u\  and  uj)  are  the 
Gaussian  ordinates  for  uj  and  rp  integrals,  respectively.  Implicit  in  Eq.  (2.36)  is  the 
implied  dependence  of  0  on  uj]  and  u\  (recall  Eq.  (2.24).  The  N-th  order  Gaussian 
ordinates  for  the  u  integral  are  [Ho75] 


WJ 


=  L^aZj  +  (1  +  aJ  (2.37) 

and  for  the  ip  integral 

"  2  "J  T  2 


Wi  =  £Zi  +  £,  (2.38) 


where  Zj   and   Zi   are  the  N  zeros  of  the  Legendre  polynomial  Pvr(Z).     The 
quadrature  weights  for  the  u  integral  are  given  by 

w.  =  (1  ~  <*>)  W.  t  (2.39) 

and  for  the  ij)  integral 

Wi  =  |Wi,  (2.40) 

where  zi  and  zj  are  evaluated  using 


2  1 1  — Z  ? 


(N+1)2   [PN+l[Zi]] 


(2.41) 
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3.    A  COMPOSITE  SKYSHINE  METHOD 

A  rigorous  treatment  of  the  shielded  silo  problem  requires  the  solution  of  a 
two-  or  three-dimensional  transport  equation.  Such  transport  solutions  tend  to  be 
computationally  expensive,  thereby  limiting  the  number  of  different  cases  which 
can  be  explored  in  design  studies.  Clearly,  there  is  a  need  for  inexpensive,  albeit 
approximate,  methods  for  estimating  the  skyshine  dose  caused  by  a  shielded 
source.  The  method  used  in  MicroSkyshine  (exponential  attenuation  with  buildup 
correction  for  radiation  penetrating  the  shield)  is  such  an  approximate  method. 
The  validity  of  the  buildup  factor  method,  however,  is  not  well  established. 
Indeed,  the  inability  of  the  buildup  factor  method  to  estimate  the  emergent  energy 
and  angular  distribution  of  photons  escaping  the  shield  should  call  for  some  caution 
when  using  the  MicroSkyshine  method.  This  skepticism  is  especially  important  for 
thick  shields  where  most  of  the  photons  leaving  the  shield  have  undergone 
interactions  in  the  shield. 

In  this  chapter,  a  composite  method  to  treat  the  shielded  skyshine  problem  is 
developed.  The  composite  method  originally  proposed  by  Keck  and  Herchenroder 
[Ke82],  first  uses  a  one-dimensional  transport  description  to  calculate  the  energy 
and  angular  distribution  of  photons  leaving  a  silo  shield.  Then  the  photons  leaving 
the  shield  are  treated  as  a  bare,  polyenergetic  anisotropic,  point  source  for  which 
the  skyshine  dose  at  a  given  detector  location  is  calculated  using  the  line-beam 
response  functions  developed  for  MicroSkyshine. 

With  this  composite  method,  the  effect  of  a  shielded  source,  in  principle,  can 
be  determined  more  accurately  than  with  the  Microskyshine  method.  The 
composite  method  can  then  be  used  to  assess  the  accuracy  of  MicroSkyshine's 
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approximate  treatment  of  a  shielded  source.  The  assessment  of  the  accuracy  of  the 
MicroSkyshine  method  is  the  central  focus  of  the  next  chapter.  In  this  chapter,  a 
rigorous  description  of  the  composite  method  along  with  a  test  of  it's  validity  is 
presented. 

3.1  Decoupling  the  Shield  and  Air  Transport  Calculations 

A  decoupling  of  the  shielded  skyshine  source  problem  can  be  performed  when 
the  number  of  photons  exiting  and  then  reentering  the  shield  after  scattering  in  the 
air  is  much  smaller  than  the  total  number  of  photons  exiting  from  the  shield. 
When  this  condition  is  satisfied,  the  calculation  of  the  energy  and  angular 
distribution  of  photons  penetrating  the  overhead  shield  becomes  independent  of  the 
subsequent  transport  of  photons  through  the  air  to  the  detector.  In  effect,  the 
source  structure  and  its  shield  have  a  negligible  effect  on  the  transport  of  the 
photons  through  the  air  once  the  photons  leave  the  source  structure.  This 
decoupling  of  a  shielded  skyshine  problem  into  two  separate  and  independent 
calculations  is  the  key  approximation  needed  for  the  composite  method. 

Before  developing  the  composite  skyshine  method  for  the  problem  of  a  point 
source  on  the  axis  of  a  cylindrical  silo  with  an  overhead  shield,  conditions  for  the 
validity  of  decoupling  the  shield  and  air  transport  problems  are  considered. 

3.1.1  Validity  of  Decoupling  for  Disk  Shields 

To  determine  the  validity  of  decoupling  a  shielded  skyshine  problem,  a 
related  problem  with  a  point  source  located  at  the  center  of  a  disk  shield  is 
considered.  For  this  problem  a  point  kernel  method,  similar  to  Kitazume's  [Ki68], 
will  be  used  to  estimate  the  number  of  photons  reflected  by  the  air  back  into  the 
shield. 
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To  determine  when  the  number  of  photons  reflecting  back  from  the  air  Into 
the  slab  shield  is  small,  a  simple  point  kernel  calculation  is  performed.  Figure  3.1 
illustrates  the  geometry  used  to  describe  this  reflection  problem.  The  reflected  flux 
of  single-scattered  photons  at  a  point  xp  on  the  shield  caused  by  a  point  source 
located  on  the  center  of  a  cylindrical  shield  emitting  photons  into  the  upper 
hemisphere  can  be  estimated  from 

dV  N  Z  e<rc(EA)  bp^p)  6       2  2  >  (3J) 

vol  rprs 

air 

where  eac(E,ds)  is  the  differential  scattering  cross  section,  rp  is  the  source-to- 
scattering  point  distance,  NZ  is  the  electron  density  of  the  air,  /jp  is  the  linear 
interaction  coefficient  at  the  source  energy,  //s  is  the  linear  interaction  coefficient  at 
the  scattered  energy,  and  9S  is  the  angle  through  which  the  photon  scatters. 
In  spherical  coordinates,  the  differential  scattering  volume  is 


dV  =  rP  drp  sin0  d0  d0  .  (3.2) 


Using  the  law  of  sines  with  Fig.  3.1  results  in  the  following  trigonometric  relations 


is/s'mcp  =  rd/sin0s  ,  (3.3) 


sTn^--sin(0+^s)'  (3"4) 
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source-detector 
axis 


dmax 


Figure  3.1.  Angular  and  geometrical  illustration  of  the  variables  with  a  point 

source  on  the  axis  of  a  cylindrical  shield  used  to  calculate  the 
re-entrant  flux  into  the  shield. 
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and, 

sm^  =  s\n\^+Os)  ■  (3-5) 

Solving  for  rs  and  rp  in  terms  of  r<i,  0S,  and  0  results  in 


_  rd  sin0  ,~a] 


and 


_    _        sin  [<t>  +  fls)  (o  7] 

Tp~    rd       sin#s ■  (li-'j 


The  derivative  drp  in  terms  of  9S  after  simplifying  and  substituting  for  sin2#s  is 


dr»=f!w  <3-s» 


Substituting  Eq.  (3.8)  into  Eq.  (3.2)  results  in  the  differential  volume  being  given 
by 


2    2 
dV  =  ^-^d0sd0d/?.  (3.9) 


Substitution  of  this  expression  for  dV  into  Eq.  (3.1)  then  gives  the  reflected  flux  at 
Xp  as 
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*(xP)  =  f*      ddf^dtf71     d6s  N  Z  e<7c(e,0s)  ^1^  (3.10) 

-7T/2  0  7T— 0 


rg  rg  e^pFp  e^5 

x  2    2 

rd  rp  rs 


or 


*(xP)  =  f    d/?/W'    d0s  N  Z  e^E,^)  %^ e_/ipFp  e^srs.         (3.11) 


Since,  the  integrand  in  Eq.  (3.11)  is  independent  of  the  angle  /?,  Eq.  (3.11)  can  be 
integrated  over  dp  to  give 

*(xP)  =  ^r^  f    d^>  f        ^NZ  e^c(E,^)  e"^prp  e-^5  .  (3.12) 

ZTd        J0  J    7T-0 

The  contribution  of  secondary  particles  produced  along  the  scattered  path  to  the 
reflected  flux  is  estimated  by  including  an  appropriate  exposure  buildup  factor 
B(E,  rs)  in  Eq.  (3.12).   Thus,  the  reflected  flux  of  photons  at  point  xP  is 


7T  7T 

*(xP)  =  |^  /    d0  /      d0s  e*c(E,0s)  B(Es,rs)  e^pFp  e^5.  (3.13) 

zrd      ^0  JTT-(f) 


The  fraction  F  of  the  photons  at  the  source  energy  being  reflected  back  to  the  slab 
surface  can  then  be  estimated  by  multiplying  Eq.  (3.13)  by  the  differential  slab 
area  dA,  integrating  overall  dA  on  the  shield  surface,  and  dividing  by  the  source 
strength  to  obtain 
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^p  ^     ^Q  •'q  rd      ^Q         ^  ^ 


x  e<7c(E^8)  B(E,rs)  e  ^e^  .  (3.14) 

Since  all  the  terms  in  the  integrand  are  independent  of  0,  Eq.  (3.14)  can  be 
integrated  over  ip  to  obtain 


F  =  ^  r'^drd  f^  f  V  e<7c(Eg0s)  B(E,rs)  e**'  e**  .        (3.15) 


Evaluation  of  Eq.  (3.15)  was  done  using  the  total  interaction  coefficient  data  from 
Hubbell  [Hu82].  The  total  interaction  coefficient  data  was  logarithmically 
interpolated  to  obtain  values  at  any  energy,  and  the  differential  scattering  cross 
sections  were  evaluated  using  the  Klein-Nishina  free-electron  model  given  by  Eq. 
(2.2). 

The  buildup  factor  B(E,rs)  was  approximated  by  an  infinite  medium 
exposure  buildup  factor  for  a  point  isotropic  source.  The  buildup  factor  used  in  the 
numerical  evaluation  of  Eq.  (3.15)  was  taken  as  the  geometric  progression  model 
proposed  by  Harima  [Ha83,  Ha86]  and  given  by  Eq.  (2.8).  The  integral  in  Eq. 
(3.15)  was  evaluated  using  Gaussian  numerical  quadrature,  [Ho75]  and  the  integral 
over  drd  was  divided  as  needed  until  the  fractional  change  in  the  integrand  was  less 
than  a  small  specified  value.  The  inner  integrals,  those  over  d(f>  and  d9s,  were 
evaluated  using  16-point  Gaussian  quadrature. 
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The  calculated  results  obtained  from  Eq.  (3.15)  for  photon  energies  of  6.129, 
1.25,  and  0.5  MeV  for  different  slab  radii  are  plotted  in  Fig.  3.2.  The  results 
shown  in  Fig.  3.2,  indicate  that  the  fraction  of  source  photons  reflected  back  to  the 
slab  increases  with  the  slab's  radius  and  decreases  with  increasing  photon  energy. 
For  the  energies  and  slab  shield  sizes  considered  in  this  study,  the  fraction  of  the 
particles  reentering  the  slab  is  very  small  with  the  largest  F  having  a  value  of 
4.2%.  It  is  therefore  concluded  that  the  reflection  of  particles  from  the  air  back  to 
the  shield  can  be  safely  ignored  for  photon  energies  between  0.5  and  6.2  MeV  when 
the  shield  radius  is  less  than  7m,  and  that  the  problem  decoupling  used  in  the 
composite  method  is  reasonable. 

3.2  The  Composite  Method  for  a  Shielded-Silo  Skyshine  Problem 

Whenever  a  skyshine  problem  can  be  decoupled  into  two  independent 
transport  calculations  (transport  through  the  source  structure  and  shielding  and 
the  transport  through  the  air),  it  is  always  better  to  take  advantage  of  this 
decoupling  rather  than  to  treat  the  skyshine  problem  as  a  single,  more  complex, 
transport  calculation. 

In  this  section,  the  composite  method  is  developed  for  the  particular  skyshine 
source  shown  in  Fig.  3.3.  In  this  skyshine  problem,  a  point  monoenergetic  source  is 
placed  on  the  axis  of  a  cylindrical  silo  with  very  thick  walls.  The  top  of  the  silo 
has  a  horizontal  slab  shield  through  which  most  of  the  radiation  that  eventually 
reaches  the  detector  must  escape. 

In  applying  the  composite  method  to  this  problem,  first  the  energy, 
directional  and  spatial  distribution  of  the  gamma  photons  leaving  the  silo  structure 
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Figure  3.2  Fraction   of  photons   from   a   point   source   on   the   axis   of   a 

cylindrical  shield  reflected  back  into  a  cylindrical  shield. 
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Figure  3.3         Photon  scattering  paths  from  a  point  source  into  a  silo  causing  a 
radiation  field  outside  the  silo. 
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are  determined.    Then  this  escaping  distribution  is  used  as  a  skyshine  source  to 
determine  the  skyshine  dose  at  locations  removed  from  the  silo. 

Both  the  silo  leakage  problem  and  the  air  transport  problems  are  in 
themselves  still  formidable  transport  problems  if  performed  rigorously.  However, 
additional  approximations  can  be  made  to  simplify  considerably  these  two 
problems.  In  the  sections  below,  approximations  are  introduced  which  allow  one  to 
use  a  simple  one-dimensional  transport  calculation  to  determine  the  silo  leakage 
and  to  use  the  MicroSkyshine  method  for  the  air-transport  phase. 

3.2.1   Leakage  from  the  Source  Silo 

The  calculation  of  the  energy,  direction,  and  position  of  photons  leaking  from 
the  source  silo  of  Fig.  3.3  is  a  difficult  task  that  requires  a  two-dimensional 
transport  calculation  (after  making  use  of  the  cylindrical  symmetry  for  this 
problem).  The  point  source  emits  (monoenergetic)  photons  isotropically  in  all 
directions,  and  these  photons  can  travel  back  and  forth  in  the  silo  scattering  off  the 
silo  walls,  floor,  shield  and  even  the  source  holder  (not  shown).  Moreover,  photons 
which  do  escape  from  the  silo  can  be  scattered  back  to  the  silo  from  outside  air 
scatters. 

The  present  silo  skyshine  problem  is  modeled  after  the  KSU  Benchmark 
Skyshine  Experiment  [C178]  in  which  the  silo  walls  were  much  thicker  than  the 
roof  shield  so  that  radial  photon  leakage  through  the  walls  was  negligible  compared 
to  that  through  the  roof  shield.  Consequently,  only  photon  leakage  through  the 
roof  shield  is  considered. 

However,  even  this  roof  leakage  component  is  complicated  by  the  in-silo 
scattering  which  can  occur.    But,  photons  which  scatter  one  or  more  times  inside 
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the  silo  before  reaching  the  roof  shield  will  have  lower  energies  than  photons  which 
reach  the  roof  directly  from  the  source  and  hence  have  less  chance  of  escaping 
through  the  shield.  Moreover,  even  those  few  in-silo  scattered  photons  which  do 
escape  will  have  even  lower  energies  and,  thus  be  preferentially  absorbed  in  the 
air-transport  phase.  Since  skyshine  doses  are  desired  only  at  distances  far  from  the 
source  silo,  one  may  ignore  the  leakage  of  source  photons  which  experience  in-silo 
scattering  prior  to  migrating  through  the  roof  shield. 

Thus  if  the  inside  walls,  floor,  and  source  equipment  are  considered  black, 
and  if  the  effect  of  photon  reflection  by  the  air  outside  the  silo  is  negligible  (as 
required  by  the  composite  method),  the  photon  leakage  calculation  can  be  modeled 
by  the  transport  of  photons  through  a  finite  slab  illuminated  on  the  bottom  by 
photons  coming  directly  from  the  point  source  and  a  vacuum  boundary  condition 
on  the  top  (i.e.,  no  incident  photons).  Using  the  geometry  of  Fig.  3.4,  the  angular 
photon  intensity  is  determined  from 

ft-V  *(r,E,ft)  +  //(r,E)  $(r,E,ft)  =  S(r,E,fl)  + 

f°°dE'  f  dfi'  ^(r,E'-*E,n'-n)  $(r',E,fi')  .        (3.16) 

In  cartesian  coordinate  the  transport  equation  is  given  by 


Q-n|f(r,E,f>)  +  fi.ty$p(?,E,ft)  +  fi-t«|§  (?,E,fi)  +  (i&E)  *(r,E,fi) 

=  Sv(r,E,ft)  +    f°°dE'    f  dH'  //s(r,E'->E,(V-*n)  $(r,E*,(V) 
(3.17) 
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X     =     t 


x    =    0 


Figure  3.4 


Roof-air  interface  coordinate  system  used  in  estimating  the 
source  condition  on  the  silo  shield's  top  surface. 
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With  the  appropriate  boundary  conditions,  Eq.  (3.17)  can  be  solved  (at  least 

-»      -» 
numerically)  for  $(r,E,ft).     From  the  solution  of  this  transport  problem,  the 

-»       -► 
leakage  distribution  of  photons  from  the  top  of  the  shield  jn(rs,E,fi)  may  be  found. 

This  distribution,  or  more  precisely,  the  angular  flow  rate  (per  unit  area)  out  of  the 

shield  surface  is  then  used  as  an  area  or  surface  source  for  the  air-transport  phase 

of  the  skyshine  calculation.   Specifically,  the  surface  source  is 

Sa(rs,E,ft)  =  n-n  *(?„E,fi)  =  jn(rfcE,fi)  (3.18) 

where  n  is  the  outward  normal  on  the  top  shield  surface  (i.e.,  unit  vector  along  the 
x-direction)  and  rs  is  a  position  on  the  top  shield  surface. 

3.2.2  Approximation  of  Leakage  by  an  Effective  Point  Source 

The  multidimensional  transport  calculation  of  the  leakage  surface  source 
from  Eq.  (3.17)  is  still  a  computationally  intensive  task.  However,  if  the  skyshine 
dose  is  to  be  calculated  at  a  distance  far  removed  from  the  source  silo,  the  spatial 
variation  of  the  escaping  photons  over  the  upper  shield  surface  is  unimportant.  In 
other  words,  all  escaping  photons  could  be  considered  as  coming  from  the  same 
point  on  the  top  of  the  shield. 

Thus,  the  skyshine  source  for  the  air-transport  phase  of  the  composite 
method  (for  detector  locations  far  from  the  silo)  may  be  taken  as  a  point,  albeit 
anisotropic,  source  with  strength 

Se«(E,ft)  =  /dA  jn(?«,E,fi)  ,  (3.19) 

where  the  dA  integration  is  performed  over  that  portion  of  the  shield  surface  from 
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which  photons  escape.  For  an  infinite  slab  shield  this  gives  an  infinite  integration 
range;  however,  in  practice  only  that  portion  of  the  shield  illuminated  directly  by 
the  actual  source  contributes  significantly  to  total  leakage,  and  thus  the  surface 
integration  is  extended  only  over  the  silo  opening. 

This  effective  skyshine  point  source  thus  represents  the  total  (surface 
integrated)  exit  current  from  the  top  surface  of  the  shield,  i.e., 

S«ff(E,n)  =  J0ut(t,E,ft)  e  f  dA  jn(?s,E,fi)  .  (3.20) 

3.2.3  Use  of  a  1-D  Model  for  Calculating  the  Effective  Skyshine  Source 

In  the  composite  method,  the  exact  point  of  emergence  of  photons  from  the 
top  shield  surface  is  unimportant.  From  Eq.  (3.20)  it  is  seen  that  only  the  total 
(integrated)  exit  current  Jout  is  needed  to  find  the  effective  point  skyshine  source. 
The  importance  of  this  feature  is  that  a  comparatively  simple  one-dimensional 

(slab  geometry)  transport  model  may  be  used  to  calculate  Jout  directly  for  a  point 

-*       -» 
source  covered  by  a  horizontal  slab  shield  without  first  having  to  find  jn(rs,E.ft) 

(which   generally   requires   a  two-dimensional   transport   calculation).      In   this 

section,  the  procedure  is  presented  for  reducing  the  two-dimensional  problem  to  a 

one-dimensional  one. 

For  a  homogeneous,  infinite  slab  of  thickness  t,  the  probability  that  a 

particle  entering  the  bottom  of  the  slab  shield  with  energy  E  and  direction  of  travel 

1)  will  emerge  with  energies  in  dE'  about  E'  and  directions  of  travel  in  dft1  about 

iV  is  denoted  by  T(t,  E-*E\  n^ft')  dE'  diV.    The  total  flow  rates  in  and  out  of  the 

slab  surfaces  caused  by  an  arbitrary  point  source  which  illuminates  a  finite  area  of 

the  bottom  of  the  slab  shield  (z  =  0)  are 
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J?*(0,E,6)  =  /  dA  ]+  (0,x,y,E,ft)  (3.21) 


jP*t(t,E',n')  =  /  dA  j.(h,x,y,E',n')  (3.22) 


Jptv 

where  jPt  are  the  angular  flow  rates  (per  unit  slab  surface  area)  into  and  out  of  the 
slab  and  the  integration  is  over  all  the  surface  area  through  which  photons  pass. 
Multiplying  T(t,E->E',  H->iV)  by  the  total  flow  rate  into  the  slab  will  give  the  total 
flow  rate  out  of  the  slab,  i.e., 

J*Jut(t,E',iV)  =  J?*(0,E,n)  T(t,E^E',<WV)  .  (3.23) 

It  is  this  total  exit  current  that  is  sought  in  order  to  define  the  effective 
skyshine  point  source  as  given  by  Eq.  (3.20).    However,  to  use  this  result  to  find 

Dt  "*     "* 

J^  , ,  one  must  first  calculate  T(t,E->E',  ft->Q').  Fortunately,  this  quantity  can  be 
obtained  by  a  simple  one-dimensional  transport  calculation. 

Consider,  the  same  infinite  slab  shield  uniformly  illuminated  on  the  bottom 
surface,  i.e.,  the  angular  flow  rate  j+(0,E,ft)  per  unit  surface  area  is  independent  of 

position  on  the  slab  surface.   For  this  case,  the  exit  angular  flow  rate  per  unit  area 

-» 
of  the  top  surface,  j"(t,E'ft'),  is  also  independent  of  the  y  and  z  coordinates. 

Moreover,  these  two  surface  flows  are  again  related  by 

j-(t,E'iV)  =  j+(0,E,ft)  T(t,E^E',iWV).  (3.24) 

Thus  T  can  be  found  by  performing  a  one-dimensional  transport  calculation  with 

-»  -♦ 

one  slab  face  uniformly  illuminated  by  j+(0,E,ft)  in  which  the  exit  flow  j"(t,E'ft')  is 

computed. 
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To  calculate  J^  .,  however,  one  can  bypass  the  calculation  of  T.  By 
comparing  Eqs.  (3.23)  and  Eqs.  (3.24)  it  is  seen  that  the  exit  flow  rate  for  the 
one-dimensional  transport  problem  becomes  J^  .  if  the  incident  flow  rate  is  simply 
taken  as 

J+(0,E,n)  =  J?*(0,E,ft).  (3.25) 


Dt  ~* 

Thus  to  find  J*'  .  (0,E,O),  perform  a  one-dimensional  transport  calculation  for  the 

-» 
angular  flux  density  $(z,E,!2)  in  slab  geometry  subject  to  the  boundary  conditions 


$(0,E,fi)  e    j  +  (0,E,Q)/n-Q  =  J?*(0,E,ft)/n-ft  ,   n-ft  >  0 

(3.26) 
$(t,E,0  )  =0    ,     nfl  <    0 


Then  the  calculated  flow  rate  out  of  the  top  surface,  i.e.,  j"(t,E,Q)  =  n-fi  $(t.E,Q), 

"     ~*  Dt 

n-0  >  0,  is  the  desired  total  flow  rate  JJJ  .  and,  thus,  the  effective  point  skyshine 
source  of  Eq.  (3.20)  given  by 


S^f(E,Q)  =  j-(t,E,fi)  =  n-ft  $(t,E,ft)  ,  n-fi  >  0.  (3.27) 

pi 


3.2.4  1-D  Boundary  Condition  for  a  Point  Collimated  Source 

Figure  3.3  shows  a  cylindrical  silo  with  a  cylindrical  concrete  shield  placed 
over  a  point  source  located  on  the  center  line  of  the  cylinder.  A  point  source  emits 
particles  isotropically,  so  that  the  energy-angular  dependence  of  photons  leaving 
the  source  is 

S(E,fi)  =  §fP  .  (3.2S) 
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It  is  also  assumed  that  there  are  negligible  interactions  in  the  air  as  photons  travel 
from  the  source  to  a  point  (x,y)  on  the  bottom  of  the  slab  shield.  Thus  the  angular 
flux  density  at  the  bottom  of  the  shield  (z  =  0)  is 

*(x,y,0,ft)  =  J§7  S(&  -  ft  W))  ,  (3-29) 

-» 
where  12'  is  the  unit  vector  in  the  direction  of  the  ray  from  the  source  to  the 

position  (x,y)  on  the  bottom  surface,  and  d  is  the  distance  from  the  point  source  to 

the  point  (x,y).  Integration  over  the  slab  surface  then  gives 

$tot(°'")  =  /  J§z  *& "  "')  dA  "  /dA  *(x^0^)  •  (3-3°) 

A  explicit  form  of  Eq.  (3.30)  can  be  obtained  by  using  a  polar  coordinate  system 
(origin  at  the  source)  rather  than  the  Cartesian  (x,y)  system.  Let  h  be  the 
perpendicular  distance  from  the  source  to  the  slab  (polar  axis),  (6,ip)  be  the  polar 
and  azimuthal  angles  to  the  point  (x,y),  and  r  be  the  perpendicular  distance  from 
the  polar  axis  to  the  point  (x,y).   For  this  coordinate  system  one  has 

d  =  h/cos0,  (3.31) 

dA  =  r  dr  #  ,  (3.32) 

r  =  h  tan#  ,  (3.33) 

and 

dr  =  hsec20d0.  (3.34) 
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Substitution  of  these  relations  into  Eq.  (3.30)  yields 

<fr(0,fi)  =  f     r  dr  f      j§g»  cos2^'  «5(ft-0')  d^  .  (3.35) 

J0  ^0      47rn 

Because  the  point  source  is  on  the  axis  of  the  cylinder,  the  problem  is  azimuthally 
symmetric  and  Eq.  (3.35)  reduces  to 


$(O,cos0)  =  J     Ig^r  cos20'  6{cos9-  cosfl")  dr  .  (3.36) 


Let  {J  =  cos#'  and  u  =  cos#  and  substitute  into  Eq.  (3.36)  to  get 

R 

$(0,w)  =  /     |g2  w'2  <5(a^u/)  r  dr  .  (3.37) 

Next,  substitute  r  =  h  tan#'  and  dr  =  h  sec20'  d#'  into  Eq.  (3.37)  and  simplify  to 
obtain 


$(0,w)  =  /  °5|f  w'2  £(w-w')  h  tanfl"  J^  d0'  .  (3.38) 


Simplification  and  use  of -do;'  =  sin#'  d#'  reduces  Eq.  (3.38)  to 


$(0,w)  =  /    |e  ^  ^  da;1  (3.39) 


or  after  integration 
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$(0,u;)  = 


IS  '  l  -  u  y-  "o  ' 


|     0    ,« 


(3.40) 


<  u>o 


To  relate  the  total  angular  flux  (i.e.  the  angular  flux  density  integrated  over  the 
finite  slab  area)  to  the  total  flow  rate  into  the  slab  recall  that 


n.fi$in(0,ft)  =  Jjn(0,ft) 


(3.41) 


With  azimuthal  symmetry  Eq.  (3.41)  reduces  to 


Jtn(0,O>)  =  U>  $in(0,w) 


(3.42) 


The  one-dimensional  transport  boundary  condition  in  terms  of  the  total  angular 
flow  rate  is  found  by  substituting  Eq.  (3.40)  into  Eq.  (3.42)  to  obtain 


'  S 

7^    ,   1    >    U)  >  LJo  > 


Jjn(o^)  = 


(3.43) 


0   ,  u>  <  lj0 


3.3  Transport  Calculation  of  the  Effective  Point  Skyshine  Source 

In  this  section,  the  calculational  procedure  used  to  find  the  emergent  photon 
energy  and  angular  distribution  on  the  shield's  top  surface  will  be  examined. 
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3.3.1  Use  of  KSLAB  for  Shield  Penetration  Calculations 

The  one-dimensional,  azimuthally  symmetric,  time-independent  transport. 
equation  for  an  inhomogeneous  medium  can  be  written  as  [Ry79] 


u;H  (z,E,w)  +  MZ'E)  *(Z>E^)  =  Q(^E,a;) 


+  /    dE'  /       du/ //s(z,E'-E,u/->u;)  $(z,EV)  (3.44) 

•'O  -1 


where  $(z,E,w)  is  the  angular  flux  density  (integrated  over  azimuth)  and  the  flux 
independent  source  is  Q(x,E,u;).  The  total  macroscopic  cross  section  is  /*(z.E). 
The  azimuthally  averaged  macroscopic  scattering  cross  section  is  defined  such  that 
//s(z,E'-»E,u/-»u;)  dE'  do;  is  the  probability  a  photon  of  energy  E'  and  direction  of 
travel  (J  will  scatter  into  energies  dE  about  E,  and  directions  of  travel  du  about  uj. 
Following  steps  outlined  by  Ryman  [Ry79],  this  one-dimensional  transport 
equation  can  be  reduced  with  the  multigroup  and  finite-difference  approximations 
to  a  form  suitable  for  a  numerical  solution.  The  finite  difference  form  of  the 
multi-group  discrete  ordinates  transport  equation  is 

=  Qg(zk+i>  ^ +  s§(zk+i'  ^'  t3-45' 

where  the  subscripts  have  the  ranges 

k  =  l,2,...,k,        i  =  l,2,...,N,        g=  1.2,...,  G, 
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and  where  $g(z,  ,u;i)  are  the  angular  flux  densities  in  energy  group  g  and  Qg(z,  x, 
Wi)  is  the  flux-independent  source  in  energy  group  g.  The  modal  positions  x, 
inside  the  slab  are  shown  in  Fig.  3.5,  and  the  position  x,     x  is  defined  by 


\H  =  (Z"  V^  •  <3-46> 


and  the  internode  spacing  by 

AkH  =  Vl  "  zk  •  <3-47) 

The  cell-centered  angular  flux  is  calculated  from  the  cell-edge  values  using  [Ry79] 


^g(zk,^i)  +  ^g(zk+1,Wi) 

$g(V'wi) =  — r-^ —  •  (3-48) 


The  final  term  in  Eq.  (3.45)  is  the  scatter  source  and  is  given  by  [Ry79] 

g     N 

Sg(xk+i'^  =    I    I  wj /^(^r^)  $g'(z^j)  •  (;3-49) 

g'=ij=i 

The  sum  over  energy  proceeds  from  the  highest  energy  group  (1)  to  the  energy 
group  g.  The  inner  sum  calculates  the  scatter  source  into  energy  group  g  and 
angular  direction  u\  from  energy  group  g'  and  directions  u-}. 

The  quadrature  set  which  defines  the  angular  directions  can  be  broken  into 
several  subregions  symmetric  about  the  angular  direction  u  =  0.  In  Fig.  3.6  an 
example  is  shown  that  breaks  the  direction  cosines  into  a  region  composed  of  6 
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Figure  3.5  Discrete  ordinates  coordinate  system  defined  inside  a 
slab  shield  as  used  in  the  one-dimensional  transport 
code  KSLAB  [Ry79]. 
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parts.  The  points  u\  and  u>>  are  angular  break  points.  The  reason  the  direction 
cosines  are  broken  into  smaller  subregions  is  the  need  to  treat  a  source  which 
contains  a  sharp  angular  cutoff  (i.e.  a  collimated  source).  As  an  example,  a 
monoenergetic  source  incident  on  an  infinite  slab  with  source  strength  given  by 
S/u;,  for  uc  <  uj  <  1,  and  0  for  u  <  uc  would  be  modeled  using  a  cross  section  set 
with  interval  regions  -1  <  u  <  ujc,  -ojc  <  u  <  ujc,  and  ljc  <  u  <  1. 
The  angular  direction  cosines  in  each  region  are  generated  using 


~       b-a   ,      b+a  (0  -m 

<j]  =  ~2~  wj  +  -%-  ,  (3.o0) 


and  the  angular  weights  using 


~       b— a  /o  -1  \ 

wj  =-^wj  ,  (3.01) 


where  {a>j}  are  the  standard  Gaussian  ordinates  with  associated  weights  {wj},  b  is 
the  direction  cosine  for  the  top  of  the  interval,  and  a  is  the  direction  cosine  for  the 
bottom  of  the  interval.  The  only  restrictions  for  the  subinterval  regions  are  the 
avoidance  of  a  direction  cosine  at  zero  (where  the  discrete  ordinates  equations 
becomes  indeterminate)  and  the  production  of  a  sufficiently  fine  angular  mesh  so 
that  there  is  at  least  one  nonzero  exact  kernel  cross  section  ^i  ,  (z,u>j-^Ui),  g'  =  1, 
.  .  .,g,  {ujj^Ui}  for  each  possible  group-to-group  cross  section  [Ry79]. 

The  group-to-group  cross  sections  //(E  f,-»E  ,o>i-»Uj)  are  evaluated  using  an 
exact  kernel  representation  [Mi77].  To  evaluate  the  multigroup  cross  sections,  the 
angular  flux  density  is  assumed  to  be  separable  in  energy  and  direction,  [Ry79],  i.e. 
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*(z,w,E)  £  *(z,w)  M(E)  .  (3.52) 


The  exact-kernel  multigroup  cross  sections  are  evaluated  using  [Ry79] 


E  E    , 

Jf      ft  f      a  ' 
'  dE  dE'  M(E')  Mz,E'-E,  w"-o/)  ,       (3.53) 

R  E 

g+1  ^g'+l 


where 

rV 

M  ,  =  dE'  M(E')  .  (3.54) 

V+l 


The  exact  kernel  cross  sections  were  evaluated  in  this  study  for  each  allowed 
energy  group-to-energy  group  transfer  for  all  the  direction  cosines  in  the  angular 
quadrature  set  using  the  code  PHOGROUP  [Ry79].  The  energy  group  structures 
used  in  evaluating  the  cross  sections  in  this  work  were  chosen  so  that  the  skyshine 
dose  caused  by  each  energy  group  were  relatively  equal.  The  energy  group 
structure  causing  each  energy  group  to  contribute  equally  occurred  when  the 
energy  structure  was  spaced  equally.  After  selecting  the  energy  group  structure, 
the  angular  group  structure  was  selected  to  insure  angular  coverage. 

3.3.2  Effective  Skyshine  Source  for  Shielded  Silo  Problem 

The  use  of  KSLAB  [Ry79]  to  solve  the  one-dimensional  discrete  ordinates 
transport  equation  required  a  multi-group  approximation  for  the  energy 
dependence  of  the  photons  emerging  from  the  shield.  To  use  MicroSkyshine's  beam 
response  functions  with  the  energy  group  structure  from  KSLAB,  the  photon 
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energies  were  taken  as  the  midpoint  energy  of  each  energy  group  used  in  the 
KSLAB  calculations.  The  midpoint  energy  for  each  energy  group  was  calculated 
by  the  cross  section  preparation  code  PHOGROUP  [Ry79]. 

A  change  made  in  the  output  from  KSLAB  was  to  remove  the  unscattered 
source  photons  from  the  calculated  angular  flux  density  $g(z,u;i)  which  contains 
both  scattered  and  unscattered  photons.  In  this  way,  the  contribution  of  the 
unscattered  and  scattered  photons  could  be  calculated  separately.  The  unscattered 
flux  component  penetrating  the  shield  is  readily  calculated  using  exponential 
attenuation  and  the  incident  angular  flux  on  the  bottom  surface  in  direction  uj\. 
The  scattered  angular  flux  density  in  direction  u\  at  the  top  shield  surface  is  then 
calculated  as 


<^cat(t,^i)  =  *  (t,«J  -  %^1  e-PgtM  j  u.  >  o  (3.55) 

where  t  is  the  thickness  of  the  slab  and  fig  is  the  total  macroscopic  cross  section  for 

energy  group  g. 

eff 
Finally,   the  effective   point   skyshine   source   S  .  (E  ,u>)    needed   for   the 

pt      § 

composite  skyshine  method  is  obtained  from  the  KSLAB  calculated  scattered 
angular  flux  densities  that  exit  from  the  top  shield  surface  (i.e.,  from  $ff      (t,u;), 


u  >  0)  by  using  Eq.  (3.27),  namely 


SJ;tff(Eg,u;)  =  Wi  **cat(M  ,  u  >  0  (3.56) 

The  effective  point  source  for  uncollided  photons  penetrating  the  shield  is  thus 

SJ;tff(Eg,u;)  =  §4^&L  e-tel"  ,  0  <  u  <  u0  .  (3.57) 
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In  the  composite  method,  these  two  effective  skyshine  sources  are  treated 
separately  since  they  generally  have  different  ranges  of  angular  support.  The 
uncollided  photons  are  collimated  by  the  silo  walls  and  emerge  only  in  the 
directions  u0  <  u  <  1,  while  the  scattered  photons  emerge  in  all  directions. 

An  interpolation  in  the  angular  flux  density  was  required  to  link  the  air 
transport  and  slab  transport  problems  together.  The  interpolation  was  required 
because,  in  general,  the  direction  cosines  used  by  KSLAB  will  not  correspond  to 

the   Gaussian   direction   cosines   used   to   integrate   Eq.    (2.34).      To   estimate 

eff 
S  ,  (Eg,a;i)  from  $g(wj)  it  was  necessary  to  interpolate  the  emergent  angular  flux 
pt 

densities  (calculated  at  discrete  u/j  values).  A  cubic  spline  interpolation  method 
was  used  to  perform  this  interpolation. 

A  problem  arose  when  the  spline  fitting  procedure  was  done  over  the  entire 
outward  angular  range  (0  <  u  <  1).  A  spline  fit  to  the  angular  flux  densities  from 
KSLAB  for  unscattered  and  scattered  components  of  the  flux  are  shown  in  Figs. 
3.7  and  3.8,  respectively.  The  scattered  angular  flux  in  Fig.  3.7  appears  to  be 
adequately  represented  by  the  spline  fit  through  the  data  points.  The  unscattered 
flux  fit  (Fig.  3.8)  however,  contains  spurious  oscillations  especially  in  the  region 
below  the  source's  collimation  angle.  These  oscillations  are  due  to  the  sharp  cutoff 
in  the  angular  source  at  the  source's  collimation  angle. 

A  better  fit  of  the  unscattered  flux  component  is  obtained  by  breaking  the 
region  into  two  parts  with  the  breakpoint  for  the  two  regions  being  the  source 
collimation  angle.  The  spline  fitting  procedure  was  then  applied  to  the  region 
above  and  below  the  breakpoint.  The  results,  shown  in  Fig.  3.9,  no  longer  contain 
spurious  oscillations  in  the  angular  flux  density  and  thus,  more  accurately 
represent  the  angular  flux. 
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Figure  3.7. 


Emergent  angular  flux  densities  spline  fit  for  the 
eighth  energy  group  of  the  N-16  source  for  a  source 
collimation  angle  of  120°. 
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Figure  3.8.  Emergent  angular  flux  densities  for  the  source  energy 
group  spline  fit  at  intermediate  points  (i.e.  between 
the  KSLAB  values)  over  the  outward  angular 
directions  (i.e.  out  of  the  slab  face)  for  the  N-16 
source  collimated  at  120°. 
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Figure  3.9.  Emergent  angular  flux  densities  for  the  source  energy 
group  spline  fit  in  two  regions  (one  above  the  source 
collimation  angle  and  the  other  below  it)  for  the 
N-16  source  collimated  at  120°. 
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3.4  Validation  of  the  Composite  Dose  Calculation  Method 

The  composite  method  developed  above  is  similar  to  a  two-step  method  used 
by  Keck  and  Herchenroder  [Ke82]  to  obtain  results  which  they  compared  to  the 
experimental  results  from  K-State's  benchmark  skyshine  experiment.  Keck  and 
Herchenroder  used  a  hybrid  method  composed  of  an  one-dimensional  discrete 
ordinates  code  ANISN  to  estimate  the  effective  point  source  used  by 
SKYSHINE-II.  ANISN  used  18  energy  groups  and  a  P-3  Legrendre  cross  section 
representation  [Ke82].  The  hybrid  method  gave  excellent  agreement  with  the 
K-State  Benchmark  experiment  despite  the  known  tendency  of  SKYSHINE-II  to 
over  predict  the  skyshine  dose  and  the  questionable  representation  of  the  photon 
scattering  cross  sections  resulting  from  a  low-order  Legrendre  expansion. 

The  composite  method,  as  developed  in  this  chapter,  uses  a  similar  procedure 
(i.e.  one-dimensional  discrete  ordinates  code  and  an  approximate  method  to 
estimate  the  skyshine  dose)  to  calculate  the  skyshine  dose.  The  main  differences 
are  in  the  detailed  theoretical  background  used,  the  improved  representation  of  the 
photon  cross  sections  afforded  by  the  exact  kernel  representation,  and  the  use  of 
the  improved  skyshine  response  functions  developed  for  MicroSkyshine. 

3.4.1  Benchmark  Experimental  Calculations 

The  group-to-group  cross  sections  for  Co-60  gamma  rays  were  calculated 
using  the  photon  data  base  developed  by  Biggs  and  Lighthill  [Bi72,  Bg68,  Bg72]. 
The  group-to-group  cross  section  types  considered  were  the  incoherent  component 
of  the  Compton  scatter,  the  production  of  annihilation  photons,  and  the 
photoelectric  effect.  The  group-to-group  cross  sections  were  developed  for  a 
concrete  shield  of  density  2.13  g/cm3,  for  the  concrete  composition  shown  in  Table 


61 


Table  3.1.  Material  composition  of  the  concrete  used  in 
preparing  the  group-to-group  cross  sections  used  in 
K-SLAB  [Ch87J. 


Element 

Mass  Fraction 

Element 

Mass  Fraction 

H 

5.558(-03)* 

Si 

3.151(-01) 

0 

4.981(-01) 

S 

1.283(-02) 

Na 

1.710(-02) 

K 

1.924(-02) 

Mg 

2.565(-03) 

Ca 

8.294(-02) 

Al 

4.575(-02) 

Fe 

1.240(-02) 

*read  as  5.558  x  10"3 
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Table  3.2.  Angular  directions  and  angular  weights  used  in 
calculating  the  exact  kernel  cross  sections  for  the 
Co-60  benchmark  calculations.  This  quadrature  set 
is  designed  for  a  source  collimation  angle  of  150.5 
degrees. 


direction  cosines 

w 

angular  weights 

0.9958127475 

0.01048910703 

0.9800979934 

0.01966458257 

0.9595946276 

0.01966458257 

0.9438798735 

0.01048910703 

0.9165603678 

0.05868640587 

0.8236414762 

0.1235771943 

0.6788851739 

0.1602817361 

0.5154093952 

0.1602817361 

0.3706530928 

0.1235771943 

0.2777342012 

0.05868640587 

0.2259078846 

0.07072276339 

0.1273009741 

0.1131564214 

0.02869406358 

0.07072276339 

Table  3.3.  Energy  group  ranges  and  average  energies  used  in  the 
Co-60  benchmark  calculations.  The  average  energies 
were  generated  by  PHOGROUP  [Ry79]  and  were 
used  by  SKYCALC  [Ba88]. 


Group 

Energy  Group  Ranges 

Average  Group  Energies 

No. 

(MeV) 

(MeV) 

1 

1.33  -  1.17 

1.249135 

2 

1.17  -  1.00 

1.083897 

3 

1.00  -  0.89 

0.944484 

4 

0.89  -  0.78 

0.834432 

5 

0.78  -  0.67 

0.724365 

6 

0.67  -  0.57 

0.619407 

7 

0.57  -  0.46 

0.514170 

8 

0.46  -  0.35 

0.403997 

9 

0.35  -  0.24 

0.293706 

10 

0.24  -  0.15 

0.193723 

11 

0.15  -  0.05 

0.093052 
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3.1,  and  for  the  average  energy  of  the  Co-60  gamma  photons  (1.33  and  1.17  MeV). 
Eleven  energy  groups  and  26  direction  cosines  in  6  subregions  were  used  in  the 
benchmark  group-to-group  cross  section  calculations  (see  Table  3.2  and  3.3). 
Output  from  PHOGROUP  [Ry79]  gave  the  group-to-group  scattering  cross 
sections,  the  total  group  cross  sections,  and  the  average  photon  energy  in  each 
energy  group. 

With  the  calculated  cross  sections,  a  plane  one-dimensional  transport  code, 
KSLAB  [Ry79],  was  used  for  the  benchmark  shield  thicknesses  of  21  and  42.8  cm. 
The  boundary  condition  used  in  the  discrete-ordinates  transport  code  was 
determined  using  Eq.  (3.43).  The  mess  spacing  used  in  the  discrete-ordinates 
calculation  was  chosen  to  guarantee  the  convergence  of  the  angular  flux  densities. 
The  maximum  mesh  size  allowable  is  calculated  using  [Ch87] 


Axmax  =  ^fsia  (3.58) 


where  o;min  is  the  direction  cosine  closest  to  the  origin  and  //g  is  the  total  group 
cross  section  with  the  largest  value.  The  inner  iteration  scheme  used  in  KSLAB 
was  considered  converged  when  the  absolute  fractional  difference  between  the 
previous  angular  fluxes  and  the  newly  calculated  angular  fluxes  was  less  than  5  x 
10"6  ("point-wise"  convergence). 

With  the  transmitted  angular  fluxes  calculated  by  KSLAB,  the  skyshine 
doses  were  then  calculated  using  the  code  SKYCALC  written  by  the  author  (see 
Appendix  A).  The  calculated  skyshine  doses  were  divided  by  the  source  strength 
to  express  the  doses  in  units  of  rad/photon.  Also,  to  plot  the  SKYCALC  results 
against  the  benchmark  experimental  data,  the  source-to-detector  distance  d  was 
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expressed  in  units  of  mass  thickness  and  the  dose  in  units  of  normalized  exposure 
(R-m2/Sr).  The  mass  thickness  source  detector  distance  in  g/cnr2  was  calculated 
with 

3  =  d  p  ,  (3.59) 

where  p  was  the  density  of  the  air  in  g/cm3  used  by  SKYCALC  in  it's  calculation 
of  the  skyshine  doses.  The  normalized  exposure  is  calculated  using  [Ch84] 

Rx(d)  =  1.154  Rd(d)  (P/Afto  ,  (3.60) 

where  1.154  is  the  conversion  factor  between  dose  and  exposure,  d  is  the 
source-to-detector  distance  in  m,  and  AQ0  is  the  source's  solid  angle  of  collimation 
given  by 

Afl0  =  2;r(l-cos0o)  •  (3.61) 

The  normalized  composite  results  are  shown  in  Fig.  3.10  along  with  results 
obtained  with  a  10-group  DOT  calculation  [Fa87],  results  from  the  MicroSkyshine 
method  [Fa87,  Sh87],  and  the  experimental  results  from  the  K-State  benchmark 
experiment  [Fa87].  The  worst  agreement  between  the  composite  method  and  the 
benchmark  experimental  data  occurs  for  small  air  mass-thickness  distances  (i.e., 
close  to  the  source).  The  composite  method  using  11  energy  groups  calculated  the 
skyshine  dose  as  well  or  better  than,  the  more  sophisticated  DOT  procedure  using 
10  energy  groups.  The  composite  method  was  mostly  conservative  (i.e.,  over 
predictive)  in  its  estimation  of  the  skyshine  dose  while  the  10  group  DOT 
calculations  were  always  underpredictive.   As  might  be  expected,  the 
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Figure  3.10.        MicroSkyshine   results    (-  -  -),    Composite   Method 

results    ( ),    and    DOT    results    ( )    compared 

to  the  experimental  data  from  the  K-State 
benchmark  experiment  [CI 78]  for  shield  thicknesses 
of  21  and  42.8  cms. 


66 


MicroSkyshine  results  do  not  agree  with  experimental  values  as  well  as  either  the 
DOT  or  the  composite  method  skyshine  dose  calculations,  although  given  the 
simplicity  of  the  MicroSkyshine  method,  the  results  are  remarkably  good. 

The  hybrid  method  of  Keck  and  Herchenroder  [Ke82]  along  with  the 
experimental  data  and  the  composite  method  is  shown  in  Fig.  3.11.  Figure  3.11 
indicates  that  the  composite  method  and  the  hybrid  method  results  are  very 
similar  over  the  entire  range  plotted.  Both  the  hybrid  method  and  the  composite 
method  agree  closely  with  the  measured  experimental  results. 

The  fraction  of  the  total  dose  that  each  group  of  photons  leaving  the  source 
shield  contributes  to  the  total  dose  in  the  composite  method  is  shown  in  Table  3.4. 
As  expected,  the  lower  energy  groups  contribute  less  to  the  total  dose  as  the 
source-to-detector  distance  increases,  indicating  that  the  lower  energy  photons  are 
being  preferentially  attenuated.  Table  3.4  and  Fig.  3.11  indicate  that  the  greater 
errors  in  the  composite  method  occur  when  the  lower  energy  groups,  as  calculated 
at  the  shield-air  interface,  contribute  larger  portions  of  the  total  dose.  To  check 
the  adequacy  of  the  energy  group  and  angular  structures  used,  the  21-cm 
benchmark  problem  was  rerun  using  a  16-group  energy  structure  and  a  32-group 
angular  structure.  The  results  for  the  new  21-cm  benchmark  case  were  with  3%  of 
the  11-group  composite  method  results  indicating  that  the  group  structure  used  in 
the  1 1-group  energy  structure  was  adequate. 
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Figure  3.11.       Comparison    using    Koch    and    Herchenroder's    data 

[Ke82]   ( )   with  the  composite  method   results 

( )    and    the    K-State    benchmark    experimental 

results  [C178]. 
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Table  3.4.  Fraction  of  the  total  dose  each  energy  group 
contributes  to  the  skyshine  dose  for  source-to- 
detector  mass  thickness  for  the  21  cm  and  42.8 
Benchmark  cases. 


Energy 

Source  Detector 

Mass  Thickness  (g/cm2^ 

) 

Group 

3.125 

21.875 

40.625 

59.375 

78.125 

A  21  cm  slab  shield 

1 
2 
3 

1.08E-01 

2.74E-02 
6.34E-02 

1.50E-01 
3.94E-02 
8.84E-02 

1.84E-01 
5.03E-O2 
1.08E-01 

2.11E-01 
6.03E-02 
1.25E-01 

2.33E-01 
7.00E-02 
1.39E-01 

4 
5 
6 

4.48E-02 
4.84E-02 
5.44E-02 

6.16E-02 
6.62E-02 
7.14E-02 

7.41E-02 
7.88E-02 
8.22E-02 

8.34E-02 
8.76E-02 
8.84E-02 

9.05E-02 
9.31E-02 
9.11E-02 

7 
8 
9 

6.48E-02 
8.57E-02 
9.73E-02 

7.62E-02 
9.45E-02 
1.03E-01 

7.94E-02 
9.14E-02 
9.34E-02 

7.85E-02 
8.36E-02 
7.86E-02 

7.58E-02 
7.48E-02 
6.37E-02 

10 
11 

12 

1.24E-01 
1.64E-01 
1.18E-01 

1.05E-01 
9.22E-02 
5.17E-02 

8.40E-02 
5.34E-02 
2.06E-02 

6.36E-02 

3.24E-02 
7.50E-02 

4.59E-02 
2.05E-02 
2.54E-03 

B  42.8 

cm  slab  shield 

1 
2 
3 

6.23E-02 
5.97E-03 
5.21E-02 

8.68E-02 
8.88E-03 
7.45E-02 

1.08E-01 
1.16E-02 
9.26E-02 

1.25E-01 
1.41E-02 
1.07E-01 

1.38E-01 

1.66E-02 
1.19E-01 

4 
5 
6 

3.99E-02 

4.54E-02 
5.46E-02 

5.67E-02 
6.50E-O2 
7.57E-02 

7.01E-O2 
8.05E-02 
9.16E-02 

8.06E-02 
9.27E-02 
1.04E-01 

8.91EM32 
1.03E-01 

1.13E-01 

7 
8 
9 

6.83E-02 
9.38E-02 
1.08E-01 

8.57E-02 
1.11E-01 
1.24E-01 

9.51E-02 
1.16E-01 
1.21E-01 

1.00E-02 

1.14E-01 
1.11E-01 

1.04E-O2 
1.10E-01 
9.78E-02 

10 
11 
12 

1.41E-01 
1.91E-01 
1.40E-01 

1.29E-01 
1.16E-01 
6.65E-02 

1.11E-01 
7.30E-02 
2.89E-02 

9.15E-02 
4.83E-02 
1.15E-02 

7.21E-02 
3.35E-02 
4.28E-03 
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4.   ASSESSMENT  OF  THE  MICROSKYSHINE  METHOD 

The  MicroSkyshine  [Sh87]  method  for  calculating  the  skyshine  dose  from  a  shielded 
skyshine  source  had  only  the  K-State  benchmark  experimental  data  for  verification 
of  the  accuracy  of  the  shield  treatment.  Accounting  for  the  shield  above  the  source 
with  exponential  attenuation  and  an  infinite-medium,  isotropic-source,  buildup 
factor  is  an  unverified  approximation.  One  of  this  report's  objectives  is  to  assess 
the  accuracy  of  the  MicroSkyshine  method. 

In  this  chapter,  the  exponential-attenuation,  buildup-factor  approach 
MicroSkyshine  uses  when  a  slab  shield  is  placed  above  a  point  source  will  be 
investigated  for  it's  accuracy.  The  accuracy  of  the  MicroSkyshine  method  will  be 
determined  using  the  more  accurate  composite  method  as  a  benchmark.  The 
difference  between  the  two  methods  will  be  expressed  as 

Fraction  Difference  =  Rrcicro(d)  -  Rcomp(d)  ,4  ^ 

Rcomp(d) 

where  Rmicro(d)  is  the  skyshine  dose  calculated  by  MicroSkyshine  at  a 
source-to-detector  distance  d  and  RComp(d)  is  the  skyshine  dose  calculated  by  the 
composite  method  for  the  same  source-to-detector  distance. 

In  the  investigation  of  MicroSkyshine,  three  different  primary  photons 
energies  were  used.  The  photon  energies  were  the  6.129  MeV  photon  from  NT  6. 
the  two  primary  photons  from  Co-60,  and  a  0.5  MeV  photon.  In  addition  to  the 
three  photon  energies  used,  four  conical  source  collimation  angles  were  used, 
namely,  160,  120,  80,  and  40  degrees. 
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4.1  Nitrogen-16  Photon  Test  Case 

The  Nitrogen-16  exact-kernel  cross  sections  [Ry79]  for  concrete  and  iron 
were  generated  using  12  energy  groups  and  32  angular  directions  (see  Tables  4.1 
and  4.2).  The  material  composition  and  the  density  (2.13  g/cm  3)  of  concrete  was 
the  same  as  that  used  for  the  benchmark  case  (see  Table  3.1).  The  iron  cross 
sections  generated  for  the  6.129  MeV  photon  used  a  material  density  of  7.86  g/cm3. 
The  exact  kernel  cross  sections  generated  using  PHOGROUP  [Ry79]  contained 
angular  break  points  to  represent  source  collimation  angles  of  160,  120,  80,  and  40 
degrees. 

One-dimensional  solutions  of  the  transport  equation  using  KSLAB  [Ry79] 
were  run  for  various  iron  and  concrete  shield  thicknesses  between  0.01  and  6  mean 
free  path  (mfp)  thicknesses.  The  spatial  mesh  spacing,  chosen  to  guarantee  the 
convergence  of  the  numerical  solutions,  was  (Axm  =  0.25  cm).  The  point 
convergence  criteria  set  for  the  Nitrogen-16  test  cases  was  1  x  10"4  (i.e.,  maximum 
fraction  difference  between  old  and  new  angular  fluxes.) 

4.1.1  Nitrogen-16  Results  for  Concrete  Shields 

The  fraction  difference  between  the  composite  and  MicroSkyshine  results  for 
N-16  photons  for  a  concrete  shield  with  160-degree  source  collimation  angle  is 
shown  in  Fig.  4.1.  This  figure  shows  that  the  MicroSkyshine  method 
underestimates  the  skyshine  dose  at  the  detector  for  source-to-detector  distances 
less  than  500  m  and  for  shield  thicknesses  greater  than  4  mfp.  The  maximum  dose 
underestimation  by  MicroSkyshine  occurs  at  a  source  detector  distance  of  250  m 
(the  minimum  source-to-detector  distance  considered)  and  is  approximately  2.5 
times  too  low.  The  maximum  overestimation  by  the  MicroSkyshine  method 
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Table  4.1.     Energy  group  structure  and  average  group  energies  used  with  the  6.129  MeV  Q— 16 
photon      in  generating  the  exact   kernel  group— to— group  cross  sections  for  iron   and   a 
The  generated   cross  sections  will  support  source   collimation   angles  of   160,    120,   80,   and   40 
degrees. 


Energy  Group  Range 

(MeV) 


Average  Group  Energy 
(MeV) 


129 
5.5 
5.0 
4.5 
4.0 
3.5 
3.0 
2.5 
2.0 
1.5 
1.0 
0.5 


5.5 
5.0 
4.5 
4.0 
3.5 
3.0 
2.5 
2.0 
1.5 
1.0 
0.5 
0.05 


5.812401 
5.248438 
4.748177 
4.247850 
3.747436 
3.246904 
2.746203 
2.245244 
1.743843 
1.241495 
0.736987 
0.239269 


Table   4.2.      Angular  direciton   cosines   and   Gaussian   Weights   used   in 
generating   the  exact   kernel  group— to— group   cross  sections  for   the   N— 16 
photon     energy     6.129     MeV     in     iron     and     concrete.        The     direction     cosines 
Gaussian   weights  will   allow   angular  source   collimation   angles  of   160,    120, 
80,   and   40   degrees. 


and 


Angular  Direction  Cosines 


Angular  Weights 


0.99320326 
0.96984631 
0.94648936 
0.92012218 
0.85286853 
0.78561488 
0.74757249 
0.67824725 
0.58779719 
0.51847196 
0.47734079 
0.39230081 
0.28134737 
0.19630739 
0.13695200 
0.03669617S 


0.016752050 

0.026S032S0 

0.016752050 

0.048235605 

0.077176968 

0.048235605 

0.046272424 

0.086749797 

0.086749797 

0.046272424 

0.056761531 

0.10641438 

0.1064143S 

0.056761531 

0.0S6S240S9 

0.086824089 
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Figure  4.1.  Fractional  difference  between  the  MicroSkyshine  method  with  the 

composite  method  for  concrete  shields  of  various  mfp  thicknesses 
using  a  N-16  source  collimated  at  160°. 
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occurrs  at  a  shield  mfp  thickness  of  1  mfp  at  the  source-to-detector  distance  of 
1500  m  and  was  approximately  1.5  times  greater  than  the  composite  method  do 

The  fractional  difference  results  for  a  120-degree  source  collimation  angle  are 
shown  in  Fig.  4.2.  The  plotted  results  show  that  the  MicroSkyshine  method 
always  underestimates  the  skyshine  dose.  The  worst  agreement  occurs  at  the 
smaller  source-to-detector  distances  where  the  MicroSkyshine  method 
underestimates  by  a  factor  of  2.3.  When  the  shield  thicknesses  are  greater  than  1 
mfp,  the  MicroSkyshine  method  was  never  closer  than  a  factor  of  1.2  times  too  low. 

The  fractional  difference  results  for  the  80  and  40  degree  source  collimation 
angle  cases  are  shown  in  Fig.  4.3  and  in  Fig.  4.4,  respectively.  The  plotted  results 
show  that  the  MicroSkyshine  method  consistently  underestimates  the  response  at 
both  source  collimation  angles,  for  all  source  detector  distances  and  for  all  mfp 
shield  thicknesses.  The  MicroSkyshine  method's  results  fall  in  a  a  band  between 
1.7  times  and  2.5  times  too  low  for  shield  thicknesses  greater  than  1  mfp  at  source 
collimation  angle  of  80  degrees.  For  the  40  degree  source  collimation  angle  the 
MicroSkyshine  method's  results  fall  in  a  band  between  1.7  times  and  5.5  times  too 
low  for  shields  with  thicknesses  greater  than  1  mfp. 

The  effect  of  shield  thickness  on  the  detector  response  is  shown  in  Fig.  4.5. 
The  160  and  120  degree  source  collimation  angle  cases  indicate  that  the  agreement 
between  the  MicroSkyshine  method  and  the  composite  method  was  reasonably 
good  for  large  source  collimation  angles.  At  the  narrower  source  collimation 
angles,  the  MicroSkyshine  method  and  the  composite  method  disagree  (see  Fig. 
4.5).  As  shown  in  Fig.  4.5,  the  composite  method's  skyshine  dose  increases  with 
increasing  shield  mfp  thickness  out  to  approximately  1  mfp.    The  increase  in  the 
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Figure  4.2.  Fractional    difference    results    comparing    the    MicroSkyshine 

method  with  the  composite  method  for  concrete  shields  of  various 
mfp  thicknesses  using  a  N-16  source  collimated  at  120°. 
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80  Degree  Collimated  N-16  Source 
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Figure  4.3.  Fractional    difference    results    comparing    the    MicroSkyshine 

method  with  the  composite  method  for  concrete  shields  of  various 
mfp  thicknesses  using  a  N-16  source  collimated  at  80°. 


76 


CD 
CJ 

c 

CD 

c_ 

CD 


o 

4-J 

U 
CD 

C 


40  Degree  Collimated  N-1B  Source 


1.0 


0.75 


0.50 


0.25  - 


0.0 


-0.25 


-0.50 


-0.75 


-1.0 


0.0 


-x  S-D  distance  1500m 


-v  S-D  distance  1250m 


-a  S-D  distance  1000m 


O ©  S-D  distance  500m 

q a  S-D  distance  250m 


■e- 


-B- 


1.0     2.0     3.0     4.0 
shield  thickness  in  mfp 


5.0 


Figure  4.4. 


Fractional  difference  results  comparing  the  MicroSkyshine 
method  with  the  composite  method  for  concrete  shields  of  various 
mfp  thicknesses  using  a  N-16  source  collimated  at  40°. 
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Figure  4.5.  Plot  of  the  skyshine  dose  rate  at  1000  m  for  various  mfp  concrete 

shields  illuminated  by  N-16  gamma  photons.     The  solid  line 

( )  is  the  composite  method  and  the  dashed  line  ( )  is  the 

MicroSkyshine  method. 
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skyshine  dose  predicted  by  the  composite  method  caused  the  major  difference 
between  the  two  methods,  as  both  methods  have  approximately  the  same 
asymptotic  slope. 

The  increase  in  the  skyshine  dose  observed  for  the  composite  method  was 
caused  by  the  slab  shield  redistributing  the  photons  to  lower  energies  and  into  new 
directions  of  travel,  effects  the  MicroSkyshine  method  ignores.  In  particular, 
photons  are  redistributed  into  lower  energies  with  directions  of  travel  that  would 
bring  them  closer  to  the  detector.  This  change  in  photon  direction  was  the 
important  mechanism  leading  to  the  increased  skyshine  dose,  especially  for  small 
source  collimation  angles.  The  decrease  in  photon  energies  has  a  smaller  effect;  a 
change  from  6.129  MeV  to  1  MeV  was  shown  to  cause  only  a  small  change  in  the 
detector  response  functions  [Sh87]. 

The  MicroSkyshine  method  showed  no  similar  increase  in  the  skyshine  dose 
as  the  collimation  angle  decreases,  since  this  method  is  incapable  of  accounting  for 
the  redirection  of  photons  into  new  directions  of  travel  that  come  closer  to  the 
detector.  It  is  possible  that  the  MicroSkyshine  method  could  be  improved  by  using 
a  semi-empirical  correction  to  the  buildup  factor  in  the  shield  which  depends  on 
the  source's  collimation  angle. 

4.1.2.  Results  for  Nitrogen-16  Shielded  by  Iron 

To  test  the  effect  of  shield  material  on  the  calculated  skyshine  doses,  an  iron 
test  case  was  run  for  a  source  collimation  angle  of  160  degrees.  Fractional 
difference  results  between  the  two  methods  are  plotted  in  Fig.  4.6  and  show  the 
same  trends  as  were  evident  in  the  concrete  shield  for  a  source  collimation  angle  of 
160  degrees.  These  trends  were  observed  despite  the  fact  that  the  calculated 
skyshine  doses  were  different  (see  Table  4.3).     The  fact  that  the  shape  of  the 
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Figure  4.6.  Fractional    difference    results    comparing    the    MicroSkyshine 

method  with  the  composite  method  for  iron  shields  of  various 
mfp  thicknesses  using  a  N-16  source  collimated  at  160°. 
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Table  4.3  Comparison  of  the  normalized  SkyShine  dose  caculated  with  the 
composite  method  for  an  iron  and  a  concrete  shield  of  1  mfp 
thickness.  The  fraction  Difference  is  calculated  using  (iron-dose  - 
concrete-dose)  /concrete-dose. 


A  real 

Iron  SkyShine 

Concrete  SkyShine 

Fraction 

Density 

Dose 

Dose 

Difference 

(g/cm2) 

(rad/photon) 

(rad/photon) 

6.250 

2.60E-20 

3.66E-20 

-0.2896 

12.500 

1.05E-20 

1.38E-20 

-€.2391 

18.750 

5.12E-21 

6.40E-21 

-0.2000 

25.000 

2.74E-21 

3.29E-21 

-0.1672 

31.250 

1.55E-21 

1.81E-21 

-0.1436 

37.500 

9.17E-22 

1.05E-21 

-0.1267 

43.750 

5.62E-22 

6.32E-22 

-0.1108 

50.000 

3.54E-22 

3.93E-22 

-0.0992 

56.250 

2.29E-22 

2.52E-22 

-0.0913 

62.500 

1.51E-22 

1.65E-22 

-0.0848 

68.750 

1.02E-22 

1.10E-22 

-0.0727 

75.000 

6.96E-23 

7.51E-23 

-0.0732 

81.250 

4.84E-23 

5.20E-23 

-€.0692 

87.500 

3.41E-23 

3.65E-23 

-0.0658 

93.750 

2.43E-23 

2.60E-23 

-0.0654 

100.000 

1.76E-23 

1.87E-23 

-0.0588 

106.250 

1.28E-23 

1.36E-23 

-0.0588 

112.500 

9.41E-24 

9.98E-24 

-0.0571 

118.750 

6.97E-24 

7.39E-24 

-0.0568 

125.000 

5.20E-24 

5.51E-24 

-0.0563 

131.250 

3.90E-24 

4.10E-24 

-0.0488 

137.500 

3.00E-24 

3.10E-24 

-O.0323 

143.750 

2.20E-24 

2.40E-24 

-0.0833 

150.000 

1.70E-24 

1.80E-24 

-O.0556 
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fraction  difference  curves  are  the  same  despite  the  different  skyshine  doses  suggi 
that  the  fraction  difference  curves  generated  for  concrete  can  be  used  for  other 
materials. 

4.2.   Cobalt-60  Photon  Test  Cases  in  Concrete 

The  Co-60  cross  sections  for  concrete  were  generated  using  12  energy  groups 
and  28  angular  directions  (see  Tables  4.4  and  4.5).  The  concrete  composition  was 
the  same  as  that  used  in  the  benchmark  case  (see  Table  3.1).  The  concrete  cross 
sections  were  generated  with  a  density  of  2.13  g/cm3.  The  cross  sections  generated 
for  concrete  allowed  source  collimation  angles  of  160,  120,  80,  and  40  degrees. 

The  energy-angle  distribution  of  photons  emerging  from  the  concrete  shield 
was  calculated  for  the  different  collimation  angles  and  different  shield  thicknesses 
using  K-SLAB.  The  convergence  criteria  used  on  all  the  runs  with  Co-60  photons 
was  5  x  10"5.  The  input  source  for  KSLAB  used  two  energy  groups  to  represent 
separately  the  1.33  and  1.17  MeV  photons  emitted  by  Co-60  decay. 

4.2.1  Results  for  Co-60  Photons  and  Concrete  Shields 

The  fractional  differences  versus  shield  thickness  are  plotted  in  Fig.  4.7  and 
4.8  for  Co-60  calculations  using  160  and  120  degree  source  collimation  angles, 
respectively.  From  these  results  it  is  seen  that  increasing  the  source-to-detector 
distance  increases  the  amount  by  which  the  MicroSkyshine  method  overpredicts 
skyshine  doses.  Figures  4.7  and  4.8  also  show  that,  as  shield  thickness  increases 
above  1  mfp,  the  amount  of  the  overestimation  or  underestimation  by 
MicroSkyshine  decreases.  The  MicroSkyshine  error  for  a  source  collimation  angle 
of  160  degrees  ranges  from  a  2.5  overestimation  factor  to  an  underestimation  by  a 
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Table  4.4  Energy  group  sturcture  and  anverage  group  energies  used  with  the 
two  Co-60  photons  to  generate  the  exact  kernel  group-to-group 
cross  sections  for  iron  and  concrete.  The  generated  cross  sections 
will  support  source  collimation  angles  of  160,  120,  80,  and  40 
degrees. 


Energy  Group  Range 

Average  Group  Energy 

(MeV) 

(MeV) 

1.33  -  1.25 

1.289790 

1.25  -  1.00 

1.122683 

1.00  -  0.90 

0.949576 

0.90  -  0.80 

0.849537 

0.80  -  0.72 

0.759677 

0.72  -  0.63 

0.674551 

0.63  -  0.54 

0.584498 

0.54  -  0.45 

0.494427 

0.45  -  0.36 

0.404329 

0.36  -  0.27 

0.314181 

0.27  -  0.18 

0.223901 

0.18  -  0.09 

0.132712 

Table  4.5  Angular  direction  cosines  and  Gaussian  weights  used  to  generate 
the  exact  kernel  group-to-group  cross  sections  for  the  Co-60 
photon  energies  of  1.33  and  1.17  Mev  in  concrete.  The  direction 
cosines  and  Gaussian  weights  will  allow  angular  source  collimation 
angles  of  160,  120,  80,  and  40  degrees. 


Angular  Direction  Cosines 


Angular  Weights 


0.9932032580 
0.9698463105 
0.9464893630 
0.9201221820 
0.8528685320 
0.7856148820 
0.7360607943 
0.6330222215 
0.5299836517 
0.4632196062 
0.3368240889 
0.2104285716 
0.1369519990 
0.0366961778 


0.01675204978 

0.02680327964 

0.01675204978 

0.04823560492 

0.07717696787 

0.04823560492 

0.07390123423 

0.1182419748 

0.07390123423 

0.09065328401 

0.1450452544 

0.09065328401 

0.08682408885 

0.08682408885 
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Figure  4.7. 


Fractional  difference  results  comparing  the  MicroSkyshine 
method  with  the  composite  method  for  concrete  shields  of  various 
mfp  thicknesses  using  a  Co-60  point  source  collimated  at  160°. 
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120  Degree  Collimated  Co-60  Source 


CD 
U 

CD 

c_ 

CO 


c 
o 

-4-> 

(J 
CD 

C_ 


1.0 


0.75 


0.50 


0.25 


0.0 


-0.25 


-0.50 


■0.75 


1.0 


0.0 


i    i    r 


-x  S-D  distance  1200m 


v v  s-D  distance  1000m 

a a  S-D  distance  800m 

o €>  S-D  distance  600m 

G 0  S-D  distance  400m 


q a  S-D  distance  200m 

I I I I I I 


1.0     2.0     3.0     4.0 
mfp  shield  thickness 


5.0 


Figure  4.8.  Fractional  difference  results  comparing  the  MicroSkyshine 
method  with  the  composite  method  for  concrete  shields  of  various 
mfp  thicknesses  using  a  Co-60  point  source  collimated  at  120°. 
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factor  of  1.5  times.    The  fraction  error  for  the  120  degree  source  collimatiorj  angle 

ranges  from  an  overestimation  by  1.8  times  to  an  underestimation  by  1.1  times. 

The  fractional  differences  for  the  80  degree  source  collimation  angle  arc 
shown  in  Fig.  4.9.  This  figure  indicates  that  the  MicroSkyshine  method  results  and 
the  composite  method  results  agree  reasonably  well  with  each  other  (within  25%) 
at  all  source-to-detector  distances  and  at  all  mfp  shield  thicknesses. 

The  fractional  differences  for  a  40  degree  source  collimation  angle  are  shown 
in  Fig.  4.10  and  it  is  seen  that  the  MicroSkyshine  results  are  within  a  factor  of  2 
for  all  source-to-detector  distances  and  for  all  the  shield  thicknesses  investigated. 
In  general,  as  the  source-to-detector  distance  increases,  the  amount  by  which  the 
MicroSkyshine  method  underpredicts  increases.  The  worst  agreement  between  the 
MicroSkyshine  method  and  the  composite  method  at  the  40  degree  source 
collimation  angle  occurs  at  source-to-detector  distance  of  1200  m  and  the 
agreement  was  no  worse  than  a  1.8  factor  underestimation. 

The  40  and  80  degree  source  collimation  angle  cases  show  an  increase  in  the 
skyshine  dose  as  the  mfp  shield  thickness  increased.  The  detector  response  at  600 
m  for  all  four  source  collimations  is  shown  in  Fig.  4.11.  The  increase  in  the 
skyshine  dose  as  the  shield  thickness  increases  is  much  less  pronounced  for  the 
Co-60  photons  than  for  the  N-16  photons.  The  smaller  increase  in  the  skyshine 
dose  indicates  that  significant  preferential  attenuation  of  the  low  energy  photons 
that  emerge  from  the  shield  is  occurring  before  the  photons  reach  the  detector. 

The  MicroSkyshine  doses  for  source  collimation  angles  of  160  and  120  degrees 
decrease  more  slowly  than  the  composite  method  results  as  the  source-to-detector 
distance  is  increased.  The  slower  decrease  in  the  skyshine  dose  for  the  160  and  120 
degree  source  collimation  angles  is  caused  by  the  MicroSkyshine  method  using  the 
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Figure  4.9.  Fractional  differences  between  MicroSkyshine  and  the  composite 

method  results  for  concrete  shields  of  various  mfp  thicknesses 
using  a  Co-60  point  source  collimated  at  80°. 
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Figure  4.10.  Fractional  differences  between  MicroSkyshine  composite  method 
results  for  concrete  shields  of  various  mfp  thicknesses  using  a 
point  Co-60  source  collimated  at  40°. 
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Figure  4.11.        Plot  of  the  skyshine  dose  rate  at  600  m  for  various  concrete  shield 
thicknesses  illuminated  by  Co-60  gamma  photons.   The  solid  line 

( )  is  the  composite  method  and  the  dashed  line  ( )  is  the 

MicroSkyshine  method. 
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original  energy  of  the  photons  in  the  buildup  factor  approach  instead  of  the  correct 
lower  energy  as  is  used  in  the  composite  method.  The  less  pronounced  increase  in 
the  skyshine  doses  is  caused  by  the  change  in  the  energy  spectrum  of  the  photons. 
Faw  and  Shultis  showed  [Sh87]  that  below  1  MeV,  the  response  functions  became 
dependent  upon  the  photons  energy  especially  at  the  longer  source-to-detector 
distances. 

4.3.  0.5  MeV  Photon  Test  Cases  in  Concrete 

Cross  sections  generated  for  0.5  MeV  photons  used  9  energy  groups  and  20 
angular  groups.  The  energy  and  angular  group  structures  used  are  shown  in  Tables 
4.6  and  4.7.  The  material  composition  of  the  concrete  remained  the  same  as  that 
used  for  the  Benchmark  Co-€0  calculations  (see  Table  3.1).  The  0.5  MeV  cross 
sections  were  generated  for  a  concrete  density  of  2.13  g/cm3.  The  0.5  MeV  cross 
sections  were  generated  for  source  collimation  angles  of  160,  120,  80,  and  40 
degrees. 

The  0.5  MeV  cases  were  calculated  by  MicroSkyshine  and  the  composite 
method  for  shield  thicknesses  between  0.01  mfp  and  5  mfp.  The  KSLAB  criterion 
used  for  ending  the  iteration  scheme  in  the  transport  equation  was  1  x  10"4 
fractional  difference  between  the  old  angular  fluxes  and  the  new  angular  flux 
values. 

4.3.1  Results  for  0.5  MeV  Photons  in  Concrete 

The  0.5  MeV  fractional  differences  between  the  MicroSkyshine  and  the 
composite  method  results  are  shown  in  Figs.  4.12  and  4.13  for  source  collimation 
angles  of  160  and  120  degrees,  respectively.  These  figures  show  that  the  amount  by 
which  MicroSkyshine  overestimates  increases  with  increasing  source-to-detector 
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Table  4.6.  Energy  group  structure  and  average  group  energies 
used  with  the  0.5  MeV  photons  to  generate  the  exact 
kernel  group-to-group  cross  sections  for  iron  and 
concrete.  The  generated  cross  sections  will  support 
source  collimation  angles  of  160,  120,  80,  and  40 
degrees. 


Energy  Group  Range 

Average  Group  Energy 

(MeV) 

(MeV) 

0.50  -  0.45 

0.4748176 

0.45  -  0.40 

0.4248010 

0.40  -  0.35 

0.3747805 

0.35  -  0.30 

0.3247541 

0.30  -  0.25 

0.2747182 

0.25  -  0.20 

0.2246636 

0.20  -  0.15 

0.1745596 

0.15  -  0.10 

0.1242421 

0.10  -  0.05 

0.0719813 

Table  4.7.  Angular  direction  cosines  and  Gaussian  weights  used 
to  generate  the  exact  kernel  group-to-group  cross 
sections  for  the  0.5  MeV  photons  in  concrete.  The 
direction  cosines  and  Gaussian  weights  will  allow 
angular  source  collimation  angles  of  160,  120,  80,  and 
40  degrees. 


Angular  Direction  Cosines 


Angular  Weights 


0.987255551 
0.952437070 
0.902997311 
0.802739753 
0.709823967 
0.556220476 
0.431035377 
0.242612801 
0.136952868 
0.036695310 


0.030153690 
0.030153690 
0.086824089 
0.086824089 
0.133022222 
0.133022222 
0.163175911 
0.163175911 
0.086824089 
0.086824089 
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Figure  4.12.  Fractional  differences  between  MicroSkyshine  and  the 
composite  method  results  for  concrete  shields  of  various 
thicknesses  using  a  point  .5  MeV  source  collimated  at  160°. 
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Figure  4.13.  Fractional  difference  results  comparing  the  Micro  Sky  shine 
method  with  the  composite  method  for  concrete  shields  of 
various  mfp  thicknesses  using  a  .5  MeV  source  collimated  at 
1200. 
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distance.  The  effect  of  the  shield  thickness  is  more  complicated;  the  fractional 
difference  goes  through  a  local  maximum  and  then  decreases  with  increasing  shield 
thickness.  The  MicroSkyshine  method's  maximum  overestimation  is  2.46  times  the 
composite  doses  for  160-degree  source  collimation  and  is  2.25  times  greater  for 
120-degree  collimation. 

The  fractional  differences  for  80  and  40  degree  source  collimation  angles  are 
shown  in  Figs.  4.14  and  4.15,  respectively.  Fractional  difference  plots  for  80  and  10 
degree  collimation  show  that  increasing  the  source-to-detector  distance  increases 
the  amount  by  which  the  MicroSkyshine  method  will  overestimate  the  detector 
response.  In  the  80  degree  collimation  case,  the  increase  in  shield  thickness 
initially  causes  the  fractional  difference  to  increase  until,  at  various  distances 
depending  on  the  source-to-<letector  distance,  the  fractional  difference  slowly 
begins  to  decrease.  In  the  40  degree  collimation  case,  the  fractional  difference 
increases  with  increasing  source-to-detector  distance.  The  MicroSkyshine 
method's  maximum  overestimation  for  80  degree  collimation  is  1.9  times  the 
composite  result,  and,  for  the  40  degree  collimation,  the  maximum  overestimation 
was  1.7  times  the  composite  method  result. 

Plotting  the  shield  thickness  against  the  detector  response  at  300m  results  in 
Fig.  4.16.  For  the  composite  method,  the  increase  in  the  skyshine  dose  with 
increasing  shield  thicknesses  has  virtually  disappeared.  The  MicroSkyshine 
calculated  doses  for  all  source  collimation  angles  now  show  an  increase  with 
increasing  shield  thicknesses  for  thin  shields.  The  MicroSkyshine  method 
overpredicts  the  skyshine  dose  because  the  buildup  factor  approach  with  the 
original  source  energy  does  not  take  into  account  the  large  energy  variation  in  the 
response  functions  for  low-energy  photons. 
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Figure  4.14.  Fractional  differences  between  MicroSkyshine  and  the  composite 
method  results  for  concrete  shields  of  various  mfp  thicknesses 
using  a  point  .5  MeV  source  collimated  at  80°. 
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40  Degree  Collimated  0.5  MeV  Source 
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Figure  4.15.  Fractional  differences  between  MicroSkyshine  and  the  composite 
method  results  for  concrete  shields  of  various  mfp  thicknesses 
using  a  point  .5  MeV  source  collimated  at  40°. 


96 


l.E-20 


c 

o 

4-J 

o 

1 

E- 

-21 

JZ 

TD 

CO 

f_ 

CD 

1 

.E- 

-22 

C/) 

O 

TD 

CD 

C 

-t-l 

.C 

en 

1 

.E- 

-23 

.m 

en 

"O 

CD 

N 

•M 

1 — 1 

CD 

F 

1 

.E- 

-24 

f_ 

O 

c 

l.E-25 


i 1 1 1 1 1 1 1 

q q160  collimation  angle  I 


a ^80  collimation  angle 

v v40  collimation  angle 

i i i i i i L 


0.0     1.0     2.0     3.0     4.0 
mfp  shield  thickness 


5.0 


Figure  4.16.        Plot  of  the  skyshine  doses  at  300  m  for  various  thicknesses  of 
concrete  shields  illuminated  by  .5  MeV  gamma  photons.     The 

solid  line  ( )  is  for  the  composite  method  and  the  dashed  line  (- 

— )  is  for  the  MicroSkyshine  method. 
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4.4.   Composite  Skyshinc  Itcsults. 

The  composite  method  results  display  two  interesting  characteristics.  The 
first  feature  of  interest  is  the  energy  dependence  of  the  skyshine  dose  in  the 
composite  method.  The  second  result  of  interest  for  the  calculated  skyshine  dose  is 
the  variation  shown  by  the  dose  for  different  shield  thicknesses  at  different  source 
collimation  angles. 

The  primary  result  observed  from  examining  the  skyshine  dose's  variation 
with  energy  is  that  the  energy  dependence  of  the  skyshine  dose  is  different  for  two 
ranges  of  source-to-detector  distance.  The  first  range  is  for  small 
source-to-detector  distances  less  than  250  to  300  m.  In  this  range,  the  skyshine 
dose  is  relatively  insensitive  to  the  energy  of  the  primary  photon  (see  Fig.  4.17)  for 
all  the  different  shield  thicknesses  examined  for  both  concrete  and  iron  shields. 

The  second  range  is  for  source-to-detector  distances  greater  than  250  to  300 
m.  In  this  range,  the  skyshine  dose  becomes  dependent  on  the  initial  photon 
energy.  A  sample  case  (see  Fig.  4.17),  chosen  as  representative,  shows  that  the 
skyshine  dose  for  the  far  range  increases  with  increasing  source  photon  energy. 

The  second  feature  of  interest  observed  with  the  composite  method  is  how 
the  skyshine  dose  varies  with  different  shield  thicknesses.  Figures  4.18,  4.19.  and. 
4.20  are  representative  of  the  variation  in  the  skyshine  dose  for  the  N-16  source, 
the  Co-60  source,  and  the  0.5  MeV  source,  respectively.  The  three  different 
sources  show  the  same  basic  trends  in  their  response  to  different  shield  thicknesses 
and  different  source  collimation  angles.  The  skyshine  dose,  in  all  cases,  decreases 
below  that  for  the  unshielded  case  when  the  shield  thickness  is  greater  than  1  mfp. 
When  the  shield  is  less  than  1  mfp  thick,  the  skyshine  dose  can  increase  with 
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Figure  4.17.  Example  plot  showing  the  source-energy  dependence  of  the 
composite  skyshine  dose  with  source-to-detector  distance  for  the 
sources  collimated  at  160°  and  covered  by  a  1.5  mfp  shield  (with 
respect  to  each  of  the  source  energies). 
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Figure  4.18.  Example  plot  showing  the  variation  of  the  composite  skyshine 
dose  with  various  degrees  of  source  collimation  for  a  N-16  point 
source  at  a  source-to-detector  distance  of  475  m. 
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Figure  4.19.  Example  plot  showing  the  variation  of  the  skyshine  dose  with 
various  source  collimation  angles  in  the  Co-60  source  with  a 
source-to-detector  distance  of  475  m. 
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Figure  4.20.  Example  plot  showing  the  variation  of  the  composite  skyshine 
dose  with  various  source  collimation  angles  for  the  .5  MeV  source 
with  a  source-to-detector  distance  of  475  m. 
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increasing  shield  thickness.  The  skyshine  dose  was  found  to  increase  for  all  cases 
when  the  source  collimation  angle  is  80  or  40  degrees.  By  contrast  the  skyshine 
dose  was  found  to  decrease  with  increasing  shield  thickness  for  the  120  and  160 
degree  source  collimation  angles. 
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5  CONCLUSIONS 

In  this  study,  an  approximate  method  has  been  studied  for  calculating  the 
skyshine  dose  caused  by  a  point  gamma-photon  source  inside  a  shielded  silo.  The 
approximate  method,  known  as  the  composite  method,  uses  a  one-dimensional, 
discrete-ordinates  transport  equation  to  estimate  the  angular  source  leaving  the 
top  of  the  cylindrical  shielding  slab.  The  angular  source  is  then  used  with  beam 
response  functions  from  MicroSkyshine  [Sh87]  to  calculate  the  skyshine  dose  at  the 
desired  points  of  interest. 

The  objectives  in  developing  the  composite  method  were  to  achieve  better 
accuracy  than  was  achieved  using  a  conventional  exponential  attenuation  buildup 
factor  approach  to  overhead  shielded  sources  and  to  explore  the  accuracy  of  the 
conventional  technique  for  different  energy  photons,  different  shield  thicknesses, 
and  different  source  collimation  angles. 

The  composite  method  was  capable  of  very  accurately  estimating  the 
skyshine  dose  measured  in  the  K-State  benchmark  experiment.  The  only  regions 
where  caution  may  be  needed  is  for  source-to-detector  mass  thicknesses  less  than  8 
g/cm2  and  greater  than  100  g/cm2.  Near  the  silo,  photons  scattering  inside  the  silo 
are  of  importance  causing  the  composite  method,  which  neglects  in— silo  scattered 
photons,  to  underpredict  the  skyshine  dose.  At  far  distances  from  the  silo,  the 
composite  method's  slope  was  slightly  more  negative  value  than  was  the  DOT 
method's  slope.  Thus,  for  far  distances  the  composite  method  could  underestimate 
the  actual  skyshine  dose. 

The  complexity  of  the  composite  method  and  it's  limited  ability  to  deal  with 
geometries  other  than  a  source  on  the  axis  of  a  silo  will  limit  its  usefulness  in  many 
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practical  applications.  The  composite  method,  unlike  the  MicroSkyshine  method, 
requires  considerable  computer  resources  for  generating  the  photon  group-to-group 
cross  sections  and  then  for  solving  the  one-dimensional  transport  equation.  The 
calculation  of  the  skyshine  dose,  once  the  angular  source  is  known,  is  fairly  easy 
and  can  be  done  using  a  microcomputer. 

The  second  objective  of  this  study,  namely,  the  evaluation  of 
MicroSkyshine's  method  of  approximating  a  shield,  yielded  complex  results  and 
two  interesting  observations.  When  the  source  is  shielded  by  1  mfp  or  thicker 
shield,  MicroSkyshine  can  underestimate  the  skyshine  dose.  The  underestimation 
was  largest  for  the  6.129  MeV  gamma  source  with  the  source  collimated  tightly  at 
40  degrees.  In  general,  the  smaller  the  source  collimation  angle  the  more  likely  it 
is  for  MicroSkyshine  to  underestimate  the  skyshine  dose.  On  the  other  hand,  when 
gamma  ray  energies  are  lower  and  the  source  collimation  angles  wider,  the 
buildup-factor  method  used  in  MicroSkyshine  overpredicts  the  skyshine  dose. 

The  actual  amount  by  which  the  buildup  factor  method  in  MicroSkyshine 
overestimates  or  underestimates  the  skyshine  dose  is  dependent  upon  the  photon 
energy,  the  source  detector  distance,  and  the  source's  angle  of  collimation.  The 
maximum  underprediction  observed  with  the  MicroSkyshine  method  was  a  factor 
of  5.5  times  lower  than  the  composite  method  result  at  a  source-to-detector 
distance  of  1500m  for  a  N-16  source  collimated  at  40  degrees.  The  maximum 
overprediction  observed  with  the  MicroSkyshine  method  was  a  factor  of  2.5  times 
higher  than  the  composite  method  result  at  a  source-to-detector  distance  of  600m 
for  the  0.5  MeV  source  collimated  at  160  degrees. 

The  primary  photons  of  interest  when  calculating  practical  skyshine  doses 
are  the  high  energy  N-16  photons.    The  results  of  this  study  for  N-16  photons 
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show  that,  for  broadly  collimated  sources  and  source  detector  distant  ter 

than  200  m,  the  MicroSkyshine  method  is  within  a  factor  of  2  for  shield  thick] 
up  to  6  mfps. 

Two  interesting  observations  were  made  during  the  course  of  this  study. 
First,  the  skyshine  dose  near  the  source  can  increase  with  increasing  shield 
thickness  for  some  source  collimation  angles.  Second,  the  skyshine  dose  is 
relatively  insensitive  to  the  energy  of  the  source  photons  for  source-to-detector 
distances  less  than  300m. 

The  skyshine  dose,  for  all  photon  energies  investigated,  increases  with  shield 
thicknesses  up  to  approximately  1  mfp  for  the  80  and  40  degree  source-collimation 
angles.  The  overhead  shield  increases  the  skyshine  dose  by  redirecting  photons 
into  directions  toward  to  the  detector.  In  the  narrower  collimation  cases  no 
uncollided  photons  are  heading  in  these  directions  and,  as  a  result,  the 
shield-scattered  photons  increase  the  skyshine  dose.  The  increase  in  the  skyshine 
dose  with  increasing  shield  thickness  is  more  noticeable  for  N-16  photons  than  it  is 
for  Co-60  and  0.5  MeV  photons. 

The  composite  method  results  showed  that,  for  source-to-detector  distances 
less  than  250  m  to  300  m,  the  skyshine  dose  is  relatively  insensitive  to  the  primary 
photon  energy.  The  range  for  which  the  primary  photon  energy  becomes 
important  is  dependent  upon  the  shield  thickness  and  the  source  collimation  angle. 
When  the  distance  is  greater  than  250  m  to  300  m  the  skyshine  dose  becomes  very 
dependent  upon  the  initial  photon  energy.  At  a  source-to-detector  distance  of 
1000m  the  skyshine  dose  varies  by  an  order  of  magnitude  between  each  of  the  three 
source  photons  studied. 
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5.1   Suggestions  For  Further  Study 

The  most  interesting  area  for  further  research  is  the  correction  of  the 
MicroSkyshine  method  for  the  shielded  source  cases.  The  MicroSkyshine  method 
could  be  corrected  by  calculating  empirical  correction  factors  for  different  energies, 
collimation  angles,  and  shield  thicknesses  using  the  composite  method  results  and 
the  MicroSkyshine  results.  The  correction  of  the  MicroSkyshine  calculations,  by 
such  empirical  factors,  would  apply  rigorously  only  for  the  cylindrical  silo  case 
since  the  present  composite  method  is  limited  to  this  relatively  simple  geometry. 

Another  recommendation  for  further  study  is  to  extend  the  composite 
method  to  different  geometries.  This  extension  will  require  the  calculation  of  the 
total  emergent  angular  current  from  the  source  shield.  Symmetry  about  the  source 
axis  would  still  be  required  for  the  use  of  a  one-dimensional, 
azimuthally-symmetric  transport  code  when  calculating  the  effective  skyshine 
point  source  on  the  top  of  a  slab  shield.  For  more  general  gemoetry, 
multi-dimensional,  azimuthally-dependent  transport  codes  would  be  required  to 
calculate  the  energy-direction  distribution  of  photons  escaping  from  the  sourse 
shield. 

Finally,  the  number  of  lower  energy  groups  in  the  MicroSkyshine  response 
function  set  should  be  expanded,  if  composite  techniques  are  to  be  applied  to 
problems  involving  low  energy  photon  sources.  As  Faw  and  Shultis  [Sh87]  showed, 
the  normalized  skyshine  dose  varied  the  most  rapidly  in  the  lower  energy  ranges 
(energies  below  1  MeV.).  These  energies  were  shown  to  be  of  importance  in  this 
study  especially  for  small  source-to-detector  distances.  An  increased  number  of 
lower  energy  groups  may  reduce  the  errors  observed  at  the  small 
source-to-detector  distances  in  this  study. 
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Appendix  A 

Presented  in  this  section  is  the  computer  code  SKYCALC  used  to  evaluate 
the  skyshine  dose  from  a  shielded  point  source  in  a  silo.  The  theoretical 
background  for  this  computer  code  is  presented  in  Chapters  2  and  3.  This  (ode 
was  written  for  machines  using  ANSI-standard  FORTRAN-77. 

The  formatted  input  data  needed  by  the  code  is  listed  in  the  code  comment 
lines.  Most  of  the  input  data  file  is  prepared  by  KSLAB  [Ry79].  The  computer 
code  is  designed  to  be  run  on  a  microcomputer. 

The  output  data  is  stored  in  a  data  file  whose  name  is  specified  by  the  user. 
Output  data  contains  the  skyshine  dose  (total  and  for  each  energy  group)  for  i 
discrete  radial  source-to-detector  distance.    A  second  data  file  specified  by  the 
program    user    contains    the    skyshine    dose    (total)    at    each    discrete    radial 
source-detector  distance. 

The  last  page  of  this  section  contains  an  example  input  file  for  the  21 -cm 
KSU  benchmark  experiment. 
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PROGRAM  TO  CALCULATE  THE  SKYSHINE  DOSE  DUE  TO  A  POINT  SOURCE 
LOCATED  ON  TOP  OF  A  FLAT  CYLINDRICAL  ROOF  SHIELD.  THE  SOURCE 
CAN  BE  A  FUNCTION  OF  DIRECTION  COSINES  OUT  OF  THE  SHIELD  AND 
OF  ENERGY.  PROGRAM  USES  OUTPUT  FROM  KSLAB  FORTRAN  AND  THE 
EDITSLAB  SUBROUTINE  AS  A  DRIVER. 

WRITTEN  BY  MICHAEL  S.  BASSET 
AS  PART  OF  MASTER  DEGREE  PROGRAM 

LIST  OF  INPUT  VARIABLES  WHICH  WILL  BE  NEEDED  BY  THE  CODE. 

IF  USING  EDITSLAB  SUBROUTINE  SOME  OF  THIS  WILL  BE  INITIALIZED 

IN  THE  OUT  FILE  CREATED  BY  THAT  ROUTINE. 


NW 

NGRP 

NQUAD 

THS 
YD 
YS 
WO 

RHO 

BP 

NUSCAT 

XMIN 

XMAX 

XDEL 

EKSLAB(I)  = 

WKSLAB(I)  = 

FLUX(J,I)  = 

NCOMP(I)  = 

EXEN(I)   = 
CROSS(I)  = 

SOURCE(I)  = 


NDTYPE 


NUMBER  OF  ANGULAR  GROUPS  USED  IN  K-SLAB  * 

NUMBER  OF  ENERGY  GROUPS  USED  IN  K-SLAB  * 

ORDER  OF  GAUSSIAN  QUADRATURE  TO  BE  USED  * 

(EITHER  16  OR  32  POINT)  * 

THICKNESS  OF  ROOF  SHIELDING  SLAB  IN  CM.  * 

DETECTOR  LOCATION  IN  RELATION  TO  THE  SILO  TOP  * 

SOURCE  LOCATIONN  IN  RELATION  TO  THE  SILO  TOP  * 

GEOMETRICAL  COLLIMATION  ANGLE  FOR  PARTICLES  ON  * 

THE  SILO  TOP.  * 

AIR  DENSITY  IN  G/CM~3  * 

SOURCE  COLLIMATION  BREAK  POINT  * 

NUMBER  OF  SOURCE  GROUPS  IN  K-SLAB  * 

MINIMUM  SOURCE  -  DETECTOR  DISTANCE  (M)  * 

MAXIMUM  SOURCE  -  DETECTOR  DISTANCE  (M)  * 

THE  DELTA  CHANGE  IN  X  BETWEEN  DOSE  CALCULATIONS.  * 

THE  ENERGY  GROUP  STRUCTURE  FROM  PHOGROUP  FORTRAN.  * 

THE  ANGULAR  GROUP  STRUCTURE  FROM  K-SLAB  FORTRAN  * 

THE  ANGULAR  FLUX  DENSITIES  FOR  ANGULAR  GROUP  J  * 

AND  ENERGY  GROUP  I  * 

THE  NUMBER  OF  THE  ENERGY  GROUP  CONTAINING  * 

AN  UNSCATTERED  FLUX  COMPONENT.  * 

THE  ENERGY  OF  THE  I'TH  UNSCATTERED  FLUX  COMPONENT  * 

THE  CROSS  SECTION  OF  THE  I'TH  ENERGY  GROUP  * 

COMPONENT  TAKEN  FROM  K-SLAB  * 

THE  NUMBER  OF  PARITICLES  EMITTED  PER  SECOND  FROM  * 

THE  I'TH  SOURCE  GROUP.  USED  TO  NORMALIZE  * 

TO  DOSE/RAD.  * 

THE  TYPE  OF  RESPONSE  FUNCTION  USED  WHEN  * 

*  CALCULATING  THE  SKYSHINE  RESPONSE  (  1  =  ABSORBED  * 

*  DOSE  RAD/PHOTON;  2  =  EXPOSURE  (R/PHOTON  * 

*********************************************************************** 

INTEGER  GP,GPADJ 

CIIARACTER*64  FNAME,  A (3) 

CHARACTER*79  TITLE 

COMMON  /BDATA/PRIGAM(12,30,3) ,EN(17) ,X(40) ,W(40) ,ANGLE(20) ,PI 

COMMON  /KSLAB/NGRP,NW,EKSLAB(25) ,WKSLAB(25) ,FLUX(25,25) ,NCOMP(25) , 
lNUSCAT,THS,SOURCE(25),EXEN(25),CROSS(25),BP,STOTAL,NFLAG(25), 
2DEL(25) 
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COMMON  /D0SED/B(25,25),STBEAM(32),YD,YS,END0SE(25),E,CE,GP, 
CGPADJ,EFACT,Rll0,MSAVE,BEAM(2O),WST0P,Bl,NqUAD,Al0LD,W0 
SET  INITIAL  VALUES  FOR  NCOMP 
DO  1   I  =  1,25 
1  NCOMP(I)   =  0 


OPEN  FILE  UNIT  8  FOR  CONTROLLING  DATA  INPUT 


* 

*  OPEN'S  FILE  TO  READ  IN  INITIAL  DATA  AND  THE  ANGLULAR  FLUX 

*  DENSITIES  FOR  ALL  ENERGY  GROUPS  AND  ANGULAR  DIRECTIONS. 

*  INPUT  NAME  OF  THE  DATA  CONTROLLING  FILE 

WRITE  (*,900) 

900  FORMAT(*  INPUT  DATA  FILE  NAME  ') 
READ  (*,901)  FNAME 

901  FORMAT(A) 

OPEN  (8,FILE=FNAME) 

*  INPUT  THE  NAME  AND  OPEN  THE  OUTPUT  FILE 

WRITE  (*,905) 

905  FORMAT('  Input  name  of  the  Output  file  ') 
READ  (*,906)  FNAME 

906  FORMAT(A) 

OPEN ( 10, FILE=FNAME, STATUS  =  'UNKNOWN') 

*  INPUT  THE  NAME  AND  OPEN  THE  PLOT  FILE 

WRITE(*,910) 

910  FORMAT('  Input  the  Name  of  the  Plot  File  ') 
READ (*, 911)  FNAME 

911  FORMAT(A) 

0PEN(2O,FILE=FNAME, STATUS  =  'UNKNOWN') 
WRITE(*,912) 

READ(*,*)NDTYPE 

912  FORMAT('  Chose  Type  of  Dose  response  desired' ,/,20X, ' 1  =  Rad/phot 
*on',/,20X,'2  =  R/photon') 

READ  (8,999,END=10000)TITLE 

999  FORMAT(A) 

READ  (8,1000,END=10000)  NW,NGRP,NQUAD,THS 
NW  =  NW/2 

1000  FORMAT  (3I4,E14.7) 

READ  (8,1001,END=10000)  YD,YS,WO,RHO,BP,NUSCAT 

1001  FORMAT  (2F12.8,3E14.7,I3) 

READ  (8,1002,END=10000)  XMIN,XMAX,XDEL 

1002  FORMAT  (3E14.7) 

READ  (8,1003,END=10000)  (EKSLAB(I),I  =  1,NGRP) 

1003  FORMAT  (5F12.8) 

READ  (8,1004,END=10000)  (WKSLAB(I),I  =  1,NW) 

1004  FORMAT  (5E15.8) 

READ  (8,1005,END=10000)  ((FLUX(J,I) ,1  =  1,NGRP),J  =  1,NW) 

1005  FORMAT  (5E15.8) 

IF  (NUSCAT.GT.O)  THEN 

READ  (8,1006,END=10000)  (NCOMP(I) ,EXEN(I) , CROSS (I) ,SOURCE(I) . 
CI  =  1,NUSCAT) 
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1006  FORMAT  (I3,F12.8,E16.8,F12.8) 

END  IF 

*** 

***  WRITE 'S  BACK  OUT  INPUT  DATA  CONTROLERS  TO  OUTPUT  FILE  ON  UNIT  10 
*** 

WRITE  (10,1009) 

1009  FORMAT(3(/),1X,40('*'),'  INPUT  AND  DATA  PARAMETERS  FOR  SKYCALC  ' 
C40('*')) 

WRITE(10,1010)TITLE 

1010  FORMAT(//,20X,A80) 
WRITE(10,1011)NW 

1011  FORMATtAlSX^A'NW'^X,' NUMBER  OF  ANGULAR  DIRECTIONS  OUT  OF  Til 
CE  SLAB  FACE.  ') 

WRITE(10,1012)NGRP 

1012  F0RMAT(18X,I2,4X,'NGRP',3X,,NUMBER  OF  ENERGY  GROUPS  FROM  K-SLAB.') 
VRITE(10,1013)NQUAD 

1013  F0RMAT(18X,I2,4X,'NqUAD»,2X,* NUMBER  OF  GAUSSIAN  QUADRATURE  POINTS 
CUSED  TO  CALCULATE  THE  DOSE.') 

WRITE  (10,1035)NUSCAT 
1035  F0RMATM8X, 12, 4X,'NUSCAT', IX, 'NUMBER  OF  ENERGY  GROUPS  CONTAINING 
CUNSCATTERED  GAMMAS') 
WRITE(10,1014)THS 

1014  F0RMAT(12X,F8.4,4X,,THS',4X,,SLAB  THICKNESS  IN  (CM).') 
WRITE(10,1015)YD 

1015  F0RMAT(12X, F8. 4,4X,,YDI ,5X,' SOURCE  HEIGHT  BELOW  SILO  TOP  (NEGATIVE 
C  FOR  HEIGHTS  ABOVE  SILO  TOP) . ' ) 

WRITE(10,1016)YS 

1016  F0RMAT(12X,F8. 4, 4X,'YS',5X, 'DETECTOR  HEIGHT  BELOW  SILO  TOP  (NEGATI 
CVE  FOR  HEIGHTS  ABOVE  SILO  TOP)') 

WRITE(10,1017)WO 

1017  F0RMAT(6X,E14. 7, 4X, 'WO', 5X, 'IMPOSED  MINIMUM  COSINE  OF  TIIETA  OVER  V 
CIIICH  DOSE  IS  INTEGRATED.') 

WRITE(10,1018)RHO 

1018  FORMAT (6X,E14. 7, 4X, 'RHO' ,4X, ' AIR  DENSITY  IN  (G/CM~3)  ') 
WRITE(10,1019)BP 

1019  FORMAT^X^H^^X^BP'^X/COLUMATION  ANGLE  OF  THE  SOURCE  ON  THE 
CSLAB  SURFACE  BOTTOM.') 

WRITE(10,1020)XMIN 

1020  FORMAT (12X,F9. 3, 4X, 'XMIN' ,3X, 'MINIMUM  SOURCE-DETECTOR  DISTANCE  (M) 
C) 

WRITE(10,1021)XMAX 

1021  FORMAT(12X,F10. 2, 4X,'XMAX',3X, 'MAXIMUM  SOURCE-DETECTOR  DISTANCE  (M 
C)') 

WRITE(10,1022)XDEL 

1022  F0RMAT(12X,F9. 3, 4X,'XDEL',3X, 'CHANGE  IN  SOURCE-DETECTOR  DISTANCE  B 
CETWEEN  DOSE  CALCULATIONS') 

WRITE(10,1023) 

1023  F0RMAT(/,25X,'***  GROUP  AVERAGE  ENERGIES  FROM  K-SLAB"S  CROSS  SECT 
CION  PREPARATION  CODE  ***') 

WRITE(10,1024) 
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1024  F0RMAT(5(3X, 'GROUP', 5X,'AV.   ENERGY')) 
WRITE] 10, 1025) (I ,EKSLAB(I) ,1=1 ,NGRP) 

1025  F()RMAT(5(5X,I2,3X,F12.6,1X)) 
WRITE(10,1026) 

1026  F0RMAT(/,2OX,'***  ANGULAR  QUADRATURE  POINTS  USED  BY  K-SLAB  DURING 
CCALCULATION  OF  THE  ANGULAR  FLUXES  ***') 

VRITE(10,1027) 

1027  F0RMAT(5(3X, 'GROUP', 3X,'ANG.  DIRECTION1)) 
WRITE(10,1028)(I,WKSLAB(I),I=1,NW) 

1028  F0RMAT(5(5X,I2,4X,E13.6,1X)) 
WRITE(10,1029) 

1029  FORMAT (/,25X,'***  INFORMATION  ABOUT  THE  SOURCE  ENERGY  GAMMAS  ') 
WRITE  (10,1030) 

1030  F0RMAT(5X, 'ENERGY  GP  CONTAINING' ,5X, 'ENERGY  SOURCE  GAMMAS' ,5X, 'TOT 
CAL  CROSS  SECTION ',5X,' GROUP  SOURCE  STRENGTH' ,/,5X, 'UNSCATTERED  GAM 
CMAS') 

DO  5  I  =  1,NUSCAT 
5    WRITE  (10,1031)NCOMP(I),EXEN(I),CROSS(I),SOURCE(I) 

1031  FORMAT(14X,I3,18X,F10.8,13X,E14.7,13X,F10.7,5X) 


*** 


CHOOSE  ORDER  OF  GAUSSIAN  QUADRATURE  FOR  USE  IN  INTEGRATING  THE 
NORMALIZED  DOSE. 


IF  (NQUAD.EQ.32)  THEN 
DO  10  I  =  1,16 
10      X(I)  =  -X(33  -  I) 
DO  20  I  =  1,16 

20  W(I)  =W(33  -  I) 
ELSE 

***  NQUAD  =  16  POINT  QUADRATURE 
DO  21  I  =  1,8 
X(I)  =  -X(41-I) 
VI)  =  W(41-I) 
X (1+8)  =  X(32+I> 

21  W(I+8)  =  W(32+I( 
END  IF 
WRITE  (*,*)  'FINISHED  CALCULATING  QUADRATURE  POINTS' 


*** 


CALL  THE  SUBROUTINE  TO  DECOMPOSE  THE  SCATTERED  AND  UNSCATTERED 
GROUPS  INTO  SEPARATE  GROUPS  OF  ONLY  SCATTERED  OR  UNSCATTERED 
COMPONENTS. 

CALL  UNGRP 

CALL  SPLINE  FITTING  ROUTINE  WHICH  WILL  ENABLE  ALL  ANGULAR 
DIRECTIONS  FOR  A  GIVEN  ENERGY  GROUP  TO  BE  ESTIMATED. 

CALL  SPLINE 

LOOP  OVER  DESIRED  SOURCE  DETECTOR  DISTANCES. 
DO  30  XII  =  XMIN,XMAX,XDEL 
WRITE  (*,*)  'CALCULATING  S-D  DISTANCE  ' , XII 
TOTDOS  =0.0 
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CALL  DOSE(XII) 
DO  40  I  =  1,NGRP 
40    TOTDOS  =  TOTDOS  +  ENDOSE(I) 

IF  (NDTYPE.EQ.l)  THEN 
***  WRITE  OUT  DOSE  RESPONSE 

VRITE(10,1100)XH, TOTDOS 
1100  FORMAT(/,'  SOURCE  DETECTOR  DISTANCE'  ,F9. 3, 3X, '  IN  M'J'DOSE  FROM 
CALL  GAMAS*,E14.7,'  IN  RAD/PHOTON') 
ELSE 
***  WRITE  OUT  EXPOSURE  RESPONSE 

WRITE(10,1110)XH,TOTDOS*1.154 
1110  FORMAT(/,*  SOURCE  DETECTOR  DISTANCE' ,F9. 3, 3X, '  IN  M',/,'DOSE  FROM 
CALL  GAMAS',E14.7,'  IN  R/PHOTON') 
END  IF 

NLOOP  =  INT(NGRP/8) 

IF  (M0D(NGRP,8).GT.O)  NLOOP  =  NLOOP  +  1 
WRITE(10,2000) 

2000  FORMAT (48X,'***  ENERGY  GROUP  NUMBERS  ***') 
DO  50  K  =  1, NLOOP 

NEND  =  K*8 

IF  (NGRP.LE.NEND)  NEND  =  NGRP 

IF  (NDTYPE.EQ.l)  THEN 

WRITE(10,2001)(J,J=(K-1)*8+1,NEND) 

2001  F0RMAT(16X,8(5X,I4,5X)) 
WRITE(10,2010)  (END0SE(J),J=(K-1)*8+1,NEND) 

2010  F0RMAT(3X,'D0SE  IN  GROUP' ,8(2X,E12.5)) 

WRITE ( 10 , 2020)  (ENDOSE( J) /TOTDOS , J= (K-l) *8+l , NEND) 

2020  F0RMAT(2X,'FRACT  OF  TOTAL' ,8(2X,E12.5) ,/) 
ELSE 
WRITE(10,2002)(J,J=(K-1)*8+1,NEND) 

2002  F0RMAT(19X,8(5X,I4,5X)) 
WRITE(10,2011)  (END0SE(J),J=(K-1)*8+1,NEND) 

2011  F0RMAT(3X, 'EXPOSURE  IN  GROUP' ,8(2X,E12.5)) 
WRITE(10,2021)  (END0SE(J)/T0TD0S,J=(K-1)*8+1,NEND) 

2021  F0RMAT(5X, 'FRACT  OF  TOTAL' ,8(2X,E12.5) ,/) 
END  IF 

50    CONTINUE 

WRITE(10,2030) 
2030  F0RMAT(128('_')) 

WRITE(20,2040)XII, TOTDOS 
2040  F0RMAT(F12.3,',',E15.8) 
30   CONTINUE 

GOTO  150 
10000  WRITE(10,1220) 

1220  FORMAT  ('  END  OF  INPUT  ENCOUNTERED  IN  INPUT  DATA  FILE.  CHECK  ') 
150   STOP 

END 
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*** 
*** 


***********************************************^r*^*^^^^4:**************  *  * 

*  * 

*  SUBROUTINE  TO  SEPARATE  THE  SCATTERED  AND  UNSCATTERED  FLUX  * 

*  COMPONENTS.  SUBROLTINE  WILL  NEED  THE  INPUT  SOURCE  FOR  EACH 

*  DIRECTION  COSINE  AND  THE  TOTAL  CROSS  SECTION  FOR  EACH  GROUP  * 

*  WHICH  CONTAINS  AN  UNSCATTERED  COMPONENT  OF  THE  FLUX  * 

******************************************************* 

SUBROUTINE  UNGRP 
REAL  MFP,NWFLUX 

DIMENSION  NVFLUX(25,25),ESAVE(25) 

COMMON  /KSLAB/NGRP,NW,EKSLAB(25) ,VKSLAB(25) ,FLUX(25,25) ,NC0MP(25) , 
1NUSCAT,THS,S0URCE(25),EXEN(25),CR0SS(25),BP,ST0TAL,NFLAG(25), 
2DEL(25) 
IADD  =  0 
STOTAL  =0.0 

LOOP  TROUGHT  ALL  ENERGY  GROUPS  LOOKING  FOR  GROUPS  WHICH  HAVE 
SCATTERED  AND  UNSCATTERED  FLUX  COMPONENTS. 
DO  10  I  =  1,NGRP 
CHECK  TO  SEE  IF  ENERGY  GROUP  CONTAINS  AN  UNSCATTERED  FLUX  COMPONENT 
IF  (I.EQ.NCOMP(I))  THEN 
STOTAL  =   STOTAL  +  SOURCE(I) 
NFLAG(I  +  IADD)  =   1 
NFLAG  I  +  IADD  +  1)  =   0 
ESAVE(I+IADD)  =  EXEN(IADD+1) 
ESAVE I+IADD+1)  =  EKSLAB(I) 
LOOP  THROUGH  ANGULAR  DIRECTIONS  FINDING  UNSCATTERED  FLUXES  AND 
SETTING  FLUXES  IN  ANGULAR  GROUPS  LESS  THAN  THE  SOURCE  COLLIMATION 
ANGLE  TO  ZERO. 

DO  20  J  =  1,NW 

IF  (WKSLAB(J).LT.BP)  THEN 
NWFLUX(J,I+IADD)  =  0.0 
ELSE 

MFP  =  CROSS(I)*THS/WKSLAB(J) 
NWFLUX(J,I+IADD)  =   SOURCE(I)*EXP(-MFP)/WKSLAB(J) 
END  IF 

NWFLUX(J, I+IADD+1)  =   FLUX(J,I)-NWFLUX(J,I+IADD) 
20       CONTINUE 

IADD  =  IADD+1 
ELSE 

J   NFLAG (I  +  IADD)  =  0 
ESAVE(I+IADD)  =  EKSLAB(I) 
DO  30  J  =  1,NW 

NWFLUX(J,I+IADD)  =  FLUX (J, I) 
30       CONTINUE 

END  IF 
10    CONTINUE 

NGRP  =   NGRP  +  IADD 
DO  40  I  =  1,NGRP 
DO  50  J  =  1,NW 
FLUX(J,I)  =   NWFLUX(J,I) 
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*** 
*** 


50     CONTINUE 

EKSLAB(I)  =  ESAVE(I) 
40   CONTINUE 
RETURN 

END 

*********************************************************************** 

*  * 

*  SUBROUTINE  TO  CALCULATE  THE  SPLINE  FITTING  COEFFICIENTS  WHICH  WILL  * 

*  ALLOW  FOR  ESTIMATION  OF  THE  ANGULAR  CURRENT  DENSITIES  AT  ALL      * 
ANGULAR  DIRECTIONS  REQUIRED  BY  THE  SKYSHINE  CALCULATION  PROCEDURE.  * 


* 


*  B(I,J)  =  IS  THE  SPLINE  FIT  COEFFICIENTS  FOR  ANGLUAR  DIRECTIONS;  * 

*  WHERE  I  CORESPONDS  TO  AN  K-SLAB  ANGULAR  DIRECTION  AND  * 

*  J  TO  THE  J' TH  ENERGY  GROUP.  * 

*  * 

SUBROUTINE  SPLINE 

DIMENSION  XS0L(25) 

INTEGER  GP,GPADJ 

COMMON  /BDATA/PRIGAM(12,30,3) ,EN(17) ,X(40) ,W(40) ,ANGLE(20) ,PI 

COMMON  /KSLAB/NGRP,NW,EKSLAB(25) ,WKSLAB(25) ,FLUX(25,25) ,NC0MP(25) , 
1NUSCAT,THS,S0URCE(25),EXEN(25),CR0SS(25),BP,ST0TAL,NFLAG(25), 
2DEL(25) 

COMMON  /D0SED/B(25,25),STBEAM(32),YD,YS,END0SE(25),E,CE,GP, 
CGPADJ,EFACT,RHO,MSAVE,BEAM(20),WSTOP,B1,NQUAD,A10LD,WO 

COMMON  /SPLDAT/A(0:25,0:25),BB(0:25) 
FIND  THE  SHARP  CUTOFF  ANGLE 

M  =  0 

DO  10  I  =  1,NW-1 

IF  ((WKSLAB(I).LE.BP).AND.(WKSLAB(I+1).GT.BP))  M  =  I 
10   CONTINUE 


*** 


*** 


MSAVE  =  M 

CHANGE  THE  ANGULAR  FLUX  INTO  AN  NORMALIZED  ANGULAR  CURRENT 
DO  20  J  =  1,NGRP 
DO  30  I  =  1,NW 
30     FLUX(I,J)  =  FLUX(I,J)*WKSLAB(I)/(4.0*PI)/STOTAL 
20   CONTINUE 


*** 


***  CALCULATE  THE  MATRIX  ELEMENTS  A  AND  SOURCE  VECTOR  BB  FOR  DETERMININ 
THE  CUBIC  SPLINE  FIT  COEFFICIENTS. 


*** 
*** 
*** 

*** 
*** 
*** 
*** 


LOOP  OVER  ALL  ENERGY  GROUPS 

DO  40  JE  =  1,NGRP 

CHECK  TO  SEE  IF  GROUP  COMPOSED  OF  UNSCATTERED  FLUX  COMPONENTS. 
IF  UNSCATTERED  COMPONENT  IS  BEING  PROCESSED;  BREAK  UP  INTO  TWO 
REGIONS  TO  AVOID  IRREGULARITIES  WHICH  ARISE  FROM  SPLINE  FITTING 
SHARP  CUTTOFFS  IN  THE  ANGULAR  FLUX  DENSITY. 

MM  =  MSAVE 

IF  (NFLAG(JE).EQ.l)  THEN 
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*** 


FIT  IN  REGION  LESS  THAN  THE  BREAK  POINT  BP 

IF  (M.EQ.O)  THEN  MM  =  NV 

BY  =  FLUX(1,JE)-SL0PE(1,JE)*VKSLAB(1) 

IF  (BY.LT.O)  BY  =  0 

DO  50  I  =  1,MM-1 
50    DEL(I)  =  WKSLAB(I+1)  -VKSLAB(I) 
***  SET  FIRST  MATRIX  ELEMENTS  IN  A 

A(1,0)  =  0.0 

A  (1,1)  =  2.0*WKSLAB(2)/DEL(1) 

A(l,2)  =  1.0 

BB(1)  =  6*((FLUX(2,JE)-FLUX(1,JE))/(DEL(1)**2)-(FLUX(1,JE)-BY)/ 
C(DEL(1)*W(1))) 

DO  60  I  =  2,MM-1 

A(I,I)  =  A1V(I) 

A  1,1-1)  =  DEL(I-1)/DEL(I) 

A(I,I+l)  =  1.0 
GO    BB(I)  =  BVALUE(I,JE) 

A(MM,MM+1)  =  0.0 

A(MM,MM-1)  =0.0 

A (MM, MM)  =0.0 

BB(MM)  =0.0 

CALL  TDMA(MM,XSOL) 

DO  70  I  =  1,MM-1 
70    B(I,JE)  =  XSOL(I) 

B(MM,JE)  =  0.0 


*** 


FIT  FOR  ANGULAR  DIRECTIONS  GREATER  THAN  BP 


IF  (M.GT.O)  THEN 
MM  =  M  +  1 

BY  =  FLUX(MM,JE)  -  SLOPE(MM,JE)*VKSLAB(MM) 

IF  (BY.LT.O)  BY  =  0 

DO  80  I  =  MM,NV-1 
80    DEL(I)=WKSLAB(I+1)  -  WKSLAB(I) 

A(MM,MM-1)  =  0.0 

A(MM,MM)  =  2.0*VKSLAB(MM+1)/DEL(MM) 

A(MM,MM+1)  =  1.0 

BB(MM)  =  6*((FHJX(MM+1,JE)-FLUX(MM,JE))/(DEL(MM)**2)- 
C(FLUX(MM,JE)-BY)/(DEL(MM)*WKSLAB(MM)))  ' 

DO  90  I  =  MM+1,NV-1 

A(I,I)  =  AIV(I) 

A(I,I-1)  =  DEL(I-1)/DEL(I) 

A(I,I+1)  =  1.0 
90    BB(I)  =  BVALUE(I,JE) 

DO  95  I  =  1,NW-MM 

A(I,I-1)  =  A(MM+I-l,MM+I-2) 

A  1,1)  =  A(MM+I-l,MM+I-r 

A(I,I+1)  =  A(MM+I-1,MM+I' 
95    BB(I)  =  BB(MM+I-1) 

A(NW-MM,NW-MM+1)  =0.0 
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CALL  TDMA(NV-MM+1,XS0L) 
DO  100  I  =  1,NV-MM 
100   B(MM+I-1,JE)  =  XSOL(I) 
B(NW,JE)  =  0.0 
END  IF 


*** 


*** 


*** 


SPLINE  FIT  FOR  ENERGY  GROUPS  COMPOSED  OF  SCATTERED  RADIATION 


ELSE 

BY  =  FLUX(1,JE)-SL0PE(1,JE)*VKSLAB(1) 

IF  (BY.LT.O)  BY  =  0.0 

DO  110  I  =  1,NV-1 
110   DEL(I)  =  WKSLAB(I+1)  -VKSLAB(I) 

A(1,0  =  0.0 

A  1,1  =  2*WKSLAB(2)/DEL(1) 

A(l,2)  =  1.0 

BB(1)  =  6*((FLUX(2,JE)-FLUX(1,JE))/(DEL(1)**2)-(FLUX(1,JE) 
C-BY)/(DEL(1)*VKSLAB(1))) 

DO  120  I  =  2,NV-1 

A(I,I)  =  A1V(I) 

A(I,I-1)  =   DEL(I-1)/DEL(I) 

A(I,I+1)  =   1.0 
120   BB(I)  =  BVALUE(I,JE) 

A(NV,NV+1)  =  0.0 

CALL  TDMA(NW,XSOL) 

DO  130  I  =  1 , NV-1 
130   B(I,JE)  =  XSOL(I) 

B  NV,JE)  =  0.0 

END  IF 
40    CONTINUE 

RETURN 

END 
FUNCTION  TO  CALCULATE  THE  SLOPE  BETWEEN  TWO  POINTS 

REAL  FUNCTION  SLOPE(NF,JE) 

COMMON  /KSLAB/NGRP,NW,EKSLAB(25) ,VKSLAB(25) ,FLUX(25,25) ,NC0MP(25) , 
1NUSCAT,THS,S0URCE(25),EXEN(25),CR0SS(25),BP,ST0TAL,NFLAG(25), 
2DEL(25) 

SLOPE  =  (FLUX(NF,JE)-FLUX(NF+1,JE))/(VKSLAB(NF)-VKSLAB(NF+1)) 

END 
FUNCTION  TO  CALCULATE  ONE  OF  THE  MATRIX  ELEMENTS 

REAL  FUNCTION  A1V(I) 

COMMON  /KSLAB/NGRP,NV,EKSLAB(25) ,VKSLAB(25) ,FLUX(25,25) ,NC0MP(25) , 
1NUSCAT,THS,S0URCE(25) ,EXEN(25) ,CR0SS(25) ,BP,ST0TAL,NFLAG(25) , 
2DEL(25) 

A1V  =  2*(VKSLAB(I+1)-WKSLAB(I-1))/DEL(I) 

END 
FUNCTION  TO  DETERMINE  THE  SOURCE  VECTOR  BB9I) 

REAL  FUNCTION  BVALUE(I,JE) 

COMMON  /KSLAB/NGRP,NW,EKSLAB(25),WKSLAB(25),FLUX(25,25),NC0MP(25), 

1NUSCAT,THS,S0URCE(25) ,EXEN(25) ,CR0SS(25) ,BP,ST0TAL,NFLAG(25) , 
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*** 


2DEL(25) 

RVALUE  =  ((FLUX(I+1,JE)-FLUX(I,JE))/(DEL(I)**2)-(FLUX(I,JE)- 
CFLUX(I-1,JE))/(DEL(I)*DEL(I-1)))*6.0 
END 

TRI-DIAGONAL  MATRIX  ALGORYTIIM  TO  SOLVE  SIMILTANEOUS  EQS. 


10 


SUBROUTINE  TDMA(N,X) 
DIMENSION  11(0:25),  AL(0:25)  ,X(25) 
COMMON  /SPLDAT/A(0:25,0:25),BB(0:25) 
11(0)  =  0.0 
AL(O)  =  0.0 
A(1,0)  =  0.0 
DO  10  I  =  1,N  -  1 

II(I)  =  A(I,I+l)/(A(I,n-A(I,I-l)*II(I-l)) 
AL(I)  =  BB(I)-A  I,I-1)*AL(I-1))/(A(I,I  -A(I,I-1)*H(I-1)) 
X(N-l)  =  AL(N-l) 
DO  20  I  =  N-2,1,-1 
20    X(I)  =  AL(I)-II(I)*X(I+1) 
END 

*  * 

*  CALCULATES  THE  SKYSHINE  DOSE  FOR  A  DISTANCE  XI  AND  AND  AN  ANGULAR  * 

*  CURRENT  DENSITY  OF  FLUX (I, J).  * 

*  PROGRAM  USED  LINE  BEAM  RESPONSE  FUNCTIONS  AND  IS  USABLE  FOR      * 

*  ONLY  FOR  CYLINDRICAL  GEOMETRY.  * 

*  STRBEAM(I)  IS  THE  SPLINE  FIT  CALCULATED  VALUES  FOR  THE  I  'Til      * 

*  DIRECTION  IN  THE  GAUSSIAN  QUADRATURE  SET  FOR  A  GIVE  ENERGY  GROUP.  * 

*  * 

SUBROUTINE  DOSE(Xl) 

INTEGER  GP,GPADJ 

COMMON  /BDATA/PRIGAM(12,30,3) ,EN(17) ,X(40) ,W(40) ,ANGLE(20) .PI 

COMMON  /KSLAB/NGRP,NW,EKSLAB(25),WKSLAB(25),FLUX(25.25),NC0MP(25), 
1NUSCAT,TIIS,S0URCE(25),EXEN(25),CR0SS(25),BP,ST0TAL,NFLAG(25), 
2DEL(25) 

COMMON  /D0SED/B(25,25),STBEAM(32),YD,YS,END0SE(25),E,CE,GP, 
CGPADJ,EFACT,RlI0,MSAVE,BEAM(2O),VST0P,Bl,NQUAD,Al0LD,W0 

SUMD  =  0.0 
LOOP  OVER  ALL  ENERGY  GROUPS. 

DO  10  JE  =  1,NGRP 

E  =  EKSLAB(JE) 

CALL  ENBRAK 
DETERMINE  IF  GROUP  COMPOSED  OF  UNSCATTERED  RADIATION 

IF  (NFLAG(JE).EQ.l)  THEN 
FIND  DOSE  CAUSED  BY  UNSCATTERED  RADIATION  IN  REGION  TWO 

AMP  =  (l-BP)/2.0 

AMB  =  (l+BP)/2.0 

DO  20  I  =  1,NQUAD 

WW  =  AMB  +  AMP*X(I) 
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*** 


*** 


*** 


***  CALCULATE  THE  SPLINE  CURRENT  DENSITIES  FOR  ENERGY  JE,  AND  ANGULAR 

***  DIRECTION  WW 

20    STBEAM(I)  =  VALUE(JE,WW) 

VSTOP  =  BP 

Bl  =  1 

CALL  SKYD0S(SD0SE,X1) 

D0SE2  =  SDOSE 
***  FIND  DOSE  CAUSED  BY  UNSCATTERED  RADIATION  IN  REGION  ONE 

AMP  =   (BP)/2.0 

AMB  =  (BP)/2.0 

DO  30  I  =  l.NQUAD 

WW  =  AMB  +  AMP*X(I) 
30    STBEAM(I)  =  VALUE (JE,WW) 

VSTOP  =  WO 

Bl  =  BP 

CALL  SKYD0S(SD0SE,X1) 

ENDOSE(JE)  =  D0SE2  +  SDOSE 

ELSE 
***  FIND  DOSE  CAUSED  BY  SCATTERED  GAMMA  RAYS 

AMP  =  (l-V0)/2 

AMB  =  (l+W0)/2 

DO  40  I  =  1,NQUAD 

WW  =  AMB+AMP*X(I) 
40    STBEAM(I)  =  VALUE(JE,W) 

VSTOP  =  VO 

Bl  =  1.0 

CALL  SKYD0S(SD0SE,X1) 

ENDOSE(JE)  =  SDOSE 

END  IF 

10    CONTINUE 

END 
********************************************************* 

*  * 

*  VERSION  2.0  NEVGAM  RESPONSE  FUNCTIONS  * 

*  FROM  FAV  AND  SHULTIS'S  MICROSKYSHINE  CODE  * 
K-STATE  EXPERIMENTAL  STATION  REPORT  189  * 
DOSE  CONVERSION  FACTOR  * 
ENERGY  GROUP  CONTAINING  PHOTON  ENERGY  E  * 
ADJACENT  ENERGY  GROUP  FOR  INTERPLOATION  * 

INTERPOLATION  FACTOR  =  (E-EGP)/(EGPADJ-EGP)  * 

*  * 

*  SUBROUTINE  VILL  DETERMINE  BETVEEN  VHICH  GROUP  STRUCTURE  ENERGIES  * 

*  USED  IN  THE  PRIGAM  DATA  SET  THE  ENERGY  E  FALLS.  * 

%^^^^^^^^^^^^^^^^^^%^^>f:^^^%^^^^^^^^^^^^^^^^^:^4:3tc^:^4:^:^4:**^^*********^:***** 

SUBROUTINE  ENBRAK 

INTEGER  GP,GPADJ 

COMMON  /KSLAB/NGRP,NV,EKSLAB(25) ,VKSLAB(25) ,FLUX(25,25) ,NC0MP(25) , 
1NUSCAT,THS,S0URCE(25),EXEN(25),CR0SS(25),BP,ST0TAL,NFLAG(25), 
2DEL(25) 

COMMON  /D0SED/B(25,25),STBEAM(32),YD,YS,END0SE(25),E,CE,GP, 
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* 

* 

CE 

= 

* 

GP 

= 

* 

GPADJ 

= 

* 

EFACT 

= 

CGPADJ , EFACT , RHO , MSAVE , BEAM ( 20) , WSTOP , B 1 , NQUAD , A 1  OLD , WO 

CE  =  (1.3078E-11+(RHO/0.001225)**2)*E 
IF  (E.GT.l)  THEN 

IF  (E.EQ.10.0)  THEN 
GP  =  1 

ELSE 
"GP  =  10  -  INT(E) 

END  IF 
ELSE  IF  (E.GE.0.5)  THEN 

GP  =  10 
ELSE  IF  (E.GE.0.15)  THEN 

GP  =  11 
ELSE 

GP  =  12 
END  IF 

EGP  =  EBAR(GP) 
IF  (E.GT.EGP)  THEN 
GPADJ  =  GP  -  1 
IF  (GPADJ. EQ.O)  GPADJ  =  2 
ELSE 

GPADJ  =  GP  +  1 

IF  (GPADJ. EQ. 13)  GPADJ  =  12 
END  IF 

DELTE  =  EBAR(GPADJ)  -  EGP 
IF  (DELTE. EQ.O)  DELTE  =1.0 
EFACT  =  (E-EGP) /DELTE 
END 

*  FUNCTION  EBAR  USE  TO  FINE  THE  MEAN  GROUP  ENERGY  FOR  GROUP  I 

REAL  FUNCTION  EBAR(I) 
IF  (I.LT.10)  THEN 

EBAR  =  10.5-1 
ELSE  IF  (I.LT.ll)  THEN 

EBAR  =0.75 
ELSE  IF  (I.LT.12)  THEN 

EBAR  =  0.325 
ELSE 

EBAR  =0.1 
END  IF 
END 

*  FUNCTION  USING  THE  SPLINE  FIT  COEFFICIENTS  FOR  THE  JE'TII  ENERGY   * 

*  GROUP  TO  FIND  THE  FITTED  CURRENT  VALUES  AT  THE  ORDINATES        * 

*  DIRECTION  WW  * 

REAL  FUNCTION  VALUE(JE,WW) 

INTEGER  GP, GPADJ 

COMMON  /KSLAB/NGRP,NV,EKSLAB(25) ,WKSLAB(25) ,FLUX(25,25) ,NC0MP(25) , 
1NUSCAT,THS,S0URCE(25),EXEN(25),CR0SS(25),BP,ST0TAL,NFLAG(25), 
2DEL(25) 

COMMON  /D0SED/B(25,25),STBEAM(32),YD,YS,END0SE(25),E,CE,GP, 
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10 


20 


CGPADJ,EFACT,RHO, MSAVE, BEAM(20),WSTOP,B1,NQUAD,A10LD, WO 
M  =  MSAVE 

IF  (NFLAG(JE).Eq.l)  THEN 
IF  (VW.LT.BP)  THEN 
JJ  =  M  -  1 
DO  10  I  =  1,M-1 

IF  ((WW.GE.WKSLAB(I)).AND.(WW.LE.WKSLAB(I+1)))  JJ  =  I 
IF  (WW.LT.WKSLAB(l))  JJ  =  1 
SAVE  =  FITSP(JJ,JE,WW) 
ELSE 

jj  =  Ny_i 

DO  20  I  =  M+1,NW-1 

IF  ((WW.GE.WKSLAB(I)).AND.(WW.LE.WKSLAB(I+1)))  JJ  =   I 
IF  (WW.LT.WKSLAB(M+l)  JJ  =  M+l 
SAVE  =  FITSP(JJ,JE,WW 
END  IF 
ELSE 

JJ  =  NW  -  1 
DO  30  I  =  1,NV-1 
30      IF  ((WW.GE.WKSLAB(I)).AND.(WW.LE.WKSLAB(I+1)))  JJ  =  I 
IF  (WW.LT.WKSLAB(l))  JJ  =  1 
SAVE  =  FITSP(JJ,JE,WW) 
END  IF 

IF  (SAVE.LT.O)  SAVE  =  0 
VALUE  =  SAVE 
END 
***  FUNCTION  THAT  FITS  A  SPLINE  FIT  AT  DIRECTION  WW  AND  FOR  ENERGY 
***  JE 

REAL  FUNCTION  FITSP(JJ,JE,VV) 

COMMON  /KSLAB/NGRP,NV5EKSLAB(25) ,WKSLAB(25) ,FLUX(25,25) ,NC0MP(25) , 
1NUSCAT,THS,S0URCE(25),EXEN(25),CR0SS(25),BP,ST0TAL,NFLAG(25), 
2DEL(25) 

COMMON  /D0SED/B(25,25) ,STBEAM(32) ,YD,YS,END0SE(25) ,E,CE,GP, 
CGPADJ,EFACT,R1I0, MSAVE, BEAM(20),WSTOP,B1,NQUAD,A10LD, WO 
FITSP  =   B(JJ,JE)*((VKSLAB(JJ+1)  -  WW  **3/DEL(JJ)  -  DEL(JJ) ' 


C(WKSLAB(JJ+lj  -WW))/6.0  +  B(JJ+1,JE)*((WW-WKSLAB(JJ))**3/DEL(JJ) 
C-  DEL(JJ)*(WW-WKSLAB(JJ)))/6.0  +  FLUX(JJ,JE)*(WKSLAB(JJ+1)  -  WW)/ 
CDEL(JJ)  +  FLUX(JJ+1,JE)*(WW-WKSLAB(JJ))/DEL(JJ) 


END 

*  * 

*  FORTRAN  VERSION  ADAPTED  BY  MICHAEL  S.  BASSETT  FROM  * 

*  VERSION  4.0  MICRO-SKYSHINE  WRITTEN  BY  R.  E.  FAW  AND  * 

*  J.  K.  SHULTIS.  PROGRAM  ONLY  DEALS  WITH  UNSHIELDED  * 

*  SILO  GEOMETRY.  * 

*  * 

************************************************************************ 


SUBROUTINE  SKYD0S(SD0SE,X1) 
INTEGER  GP,GPADJ 


COMMON  /BDATA/PRIGAM(12,30,3) ,EN(17) ,X(40) ,W(40) ,ANGLE(20) ,PI 
COMMON  /KSLAB/NGRP,NW," 


EKSLAB(25) ,WKSLAB(25) ,FLUX(25,25) ,NC0MP(25) , 
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1NUSCAT,TIIS,S0URCE(25)  ,EXEN(25)  ,CR0SS(25)  ,BP,ST0TAL,NFLAG(25) , 
2DEL(25) 

COMMON  /D0SED/B(25,25),STBEAM(32),YD,YS,END0SE(25),E,CE,GP, 
CGPADJ,EFACT,RHO,MSAVE,BEAM(20),VSTOP,B1, NQUAD, AlOLI),VO 
COMMON  /SIMPLE/  R,D1,D2 
AA  =  SQRT(X1*X1  +  (YD-YS)**2) 
R  =  (RHO/0.001225)*AA 
Dl  =  Xl/AA 
D2  =  (YD-YS)/AA 
***  CALCULATE  THE  INTERPOLATED  BEAM  RESPONSE  FUNCTIONS 
***  BEAM(J)  IS  THE  BEAM  RESPONSE  FUNCTION  FOR  ALL  ANGULAR  DIRECTIONS 
***  AT  THE  INTERPOLATED  ENERGY  OF  E. 
DO  10  J  =  1,20 

BINTER  =  EXP(PRIGAM(GP,J,1))*R**PRIGAM(GP,J,2)*EXP(-R* 
CPRIGAM(GP,J,3)) 
BADJ  =  EXP(PRIGAM(GPADJ,J,1))*R**PRIGAM(GPADJ,J,2)*EXP(-R* 
CPRIGAM(GPADJ,J,3)) 
10    BEAM(J)=  BINTER  +  (BADJ-BINTER)*EFACT 
***  CALCULATE  DOSE  WITH  INTERPOLATED  RESPONSE  FUNCTIONS 
A10LD  =  -4.5 

CALL  GAUSS0(O.O,PI,SD0SE) 
SDOSE  =  2.0*SDOSE  *  CE 
END 

*  * 

*  NQUAD  GAUSSIAN  INTEGRATION  TO  FIND  THE  DOSE  CAUSED  * 

*  BY  ENERGY  GROUP  JE  * 

SUBROUTINE  GAUSS0(A2,B2,ANS2) 
INTEGER  GP  GP\DJ 

COMMON  /BD ATA/PRIG AM ( 12 , 30 , 3) ,EN(17) ,X(40) ,W(40) ,ANGLE(20) ,PI 
COMMON  /D0SED/B(25,25),STBEAM(32),YD,YS,END0SE(25),E,CE,GP, 
CGPADJ,EFACT,RIIO,MSAVE,BEAM(20),VSTOP,B1, NQUAD, A10LD, WO 
COMMON  /SIMPLE/  R,D1,D2 
AMB2  =  (B2-A2)/2.0 
APB2  =  (B2+A2)/2.0 
SUM2  =  0.0 
DO  10  I  =  1, NQUAD 
PSI  =  X(I)*AMB2  +  APB2 
10    SUM2  =  SUM2  +  V(I)*GAUSSI(PSI) 
ANS2  =  SUM2  *  AMB2 
END 


*** 


INNER  INTEGRATION  DONE  BY  GAUSSIAN  QUADRATURE. 

REAL  FUNCTION  GAUSSI(PSI) 
REAL  INBEAM 

DIMENSION  BB(32),CC(32) 
INTEGER  GP,GPADJ 
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COMMON  /BDATA/PRIGAM(12,30,3) ,EN(17) ,X(40) ,V(40) ,ANGLE(20) ,PI 
COMMON  /D0SED/B(25,25),STBEAM(32),YD,YS,END0SE(25),E,CE,GP, 
CGPADJ,EFACT,RHO,MSAVE,BEAM(20),A1,B1,NQUAD,A10LD,WO 
COMMON  /SIMPLE/  R,D1,D2 
SUM1  =  0.0 

IF  ((ABS(A10LD-A1)).LT. 0.0001)  THEN 
DO  10  K  =  1,NQUAD 
D  =  57.29577951*AC0S(C0S(PSI)*CC(K)-BB(K)) 
10       SUM1  =  SUM1  +  W(K)*INBEAM(D)*STBEAM(K) 
ELSE 
AMB1  =  (Bl-Al)/2.0 
APB1  =  (Bl+Al)/2.0 
DO  20  K  =  1,NQUAD 

OMEGA  =  X(K)*AMB1  +  APB1 
CC(K)  =  SQRT(1  -  0MEGA**2)  *  Dl 
BB(K)  =  0MEGA*D2 

D  =  57.29577951*AC0S(C0S(PSI)*CC(K)-BB(K)) 
20       SUM1  =  SUM1  +  W(K)*INBEAM(D)*STBEAM(K) 
END  IF 
A10LD  =  Al 
GAUSSI  =  SUM1*AMB1 
END 


*** 


ANGULAR  INTERPOLATION  OF  BEAM  RESPONSE  FUNCTIONS 

REAL  FUNCTION  INBEAM(D) 
INTEGER  GP,GPADJ 

COMMON  /BDATA/PRIGAM(12,30,3) ,EN(17) ,X(40) ,W(40) ,ANGLE(20) ,PI 
COMMON  /D0SED/B(25,25),STBEAM(32),YD,YS,END0SE(25),E,CE,GP, 
CGPADJ,EFACT,RIIO,MSAVE,BEAM(20),VSTOP,B1,NQUAD,A10LD,WO 
IF  (D.GE.100)  THEN 
JVALUE  =  20-INT((180-D)/20) 
IF  (D.LT.ANGLE(JVALUE))  THEN 

JADJ  =  JVALUE  -  1 
ELSE 

JADJ  =  JVALUE  +  1 
END  IF 

IF  (JADJ. EQ. 21)  JADJ  =  19 
ELSE  IF  (D.GE.20)  THEN 

JVALUE  =  16  -  INT((100-D)/10) 
IF  (D.LT. ANGLE? JVALUE))  THEN 

JADJ  =  JVALUE  -  1 
ELSE 

JADJ  =  JVALUE  +  1 
END  IF 
ELSE  IF  (D.GE.15)  THEN 
JVALUE  =  8 
IF  (D.LT.ANGLE(8))  THEN 

JADJ  =  7 
ELSE 
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JADJ  =  9 
END  IF 
ELSE  IF  (D.GE.10)  THEN 
JVALUE  =  7 
IF  (D.LT.ANGLE(7))  THEN 

JADJ  =  6 
ELSE 

JADJ  =  8 
END  IF 
ELSE  IF  (D.GE.7)  THEN 
JVALUE  =  6 
IF  (D.LT.ANGLE(6))  THEN 

JADJ  =  5 
ELSE 

JADJ  =  7 
END  IF 
ELSE  IF  (D.GE.5)  THEN 
JVALUE  =  5 
IF  (D.LT.ANGLE(5))  THEN 

JADJ  =  4 
ELSE 

JADJ  =  6 
END  IF 
ELSE  IF  (D.GE.3)  THEN 
JVALUE  =  4 
IF  (D.LT.ANGLE(4))  THEN 

JADJ  =  3 
ELSE 

JADJ  =   5 
END  IF 
ELSE 

J  JVALUE  =  INT(D)  +  1 
IF  (D.LT.ANGLE(JVALUE))  THEN 

JADJ  =   JVALUE  -  1 
ELSE 

J  JADJ  =  JVALUE  +  1 
END  IF 

IF  (JADJ.EQ.O)  JADJ  =  2 
END  IF 

INBEAM  =  BEAM(JVALUE)  +  (BEAM(JADJ)  -  BEAM(JVALUE))*(D  -  ANGLE( 
CJVALUE))/(ANGLE(JADJ)  -  ANGLE(JVALUE)) 
END 

'Library  Skydata  v.  4.0,  20  March  87 


* 


* COEFFICIENTS  FOR  BEAM  RESPONSE  FUNCTIONS  —  PRIGAM(12 ,20,3) 

* Fitted  to  calculations  of  single  Klein-Nishina  scattering  and 

* pair  production  in  line  beams  followed  by  infinite  medium  buildup 

* with  the  GP  approximation  for  buildup.    Faw  k   Shultis  March  87. 

BLOCK  DATA  RRADAT 
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COMMON  /BDATA/PRIGAM(12,30 
9.5  MeV 
'DATA  ((PRIGAM(1,J,I),I  =  1 
C-8. 91568, -0.99301,0. 00248,- 
C-10. 78282, -0.96765, 0.00243 
C-12. 08827, -0.93680,0. 00252 
C-13. 56202, -0.89569, 0.00288 
C-15. 31173, -0.82844, 0.00404 
C-16. 93617, -0.71375, 0.00636 
C-18. 13527,-0.56366,0.00859 
C-18. 93738, -0.46223, 0.01000 
C-18. 76997, -0.59498, 0.01017 
C-18. 68741, -0.69790, 0.01011 
8.5  MeV 

DATA  ((PRIGAM(2,J,I),I  =  1 
C-8. 83156, -0.99313, 0.00258, 
C-10. 69426, -0.96639, 0.00253 
C-ll. 98481, -0.93399, 0.00261 
C-13. 43509, -0.89185, 0.00297 
C-15. 19481,-0.82114,0.00411 
C-16. 85502, -0.70466, 0.00639 
C-18. 03689, -0.56408,0. 00858 
C-18. 80175, -0.47461, 0.00999 
C-18. 64283, -0.60839, 0.01017 
C-18. 55371, -0.71538, 0.01010 
7.5  MeV 

DATA  ((PRIGAM(3,J,I),I  =  1 
C-8. 73586, -0.99324, 0.00271, - 
C-10. 59345, -0.96502, 0.00265 
C-ll. 86703, -0.93148, 0.00273 
C-13. 29309, -0.88770, 0.00308 
C-15. 05540, -0.81451, 0.00420 
C-16. 74823, -0.69705,0. 00643 
C-17. 91529, -0.56660, 0.00858 
C-18. 64464, -0.48929, 0.00998 
C-18. 49817, -0.62380, 0.01017 
C-18. 40158, -0.73612,0.01009 
6.5  MeV 

DATA  ((PRIGAM(4,J,I),I  =  1 
C-8 . 63775 , -0 . 99052 , 0 . 00289 ,  ■ 
C-10. 51728, -0.95374, 0.00283 
C-ll. 79500, -0.91359, 0.00292 
C-13. 20231, -0.86576, 0.00326 
C-14. 96438, -0.78976, 0.00435 
C-16. 68261, -0.67249, 0.00653 
C-17. 85307, -0.54839, 0.00864 
C-18 . 65002 ,-0 . 45863 , 0 . 01009 
C-18. 40481, -0.62280, 0.01023 
C-18 .31056 ,-0 . 73909 , 0 . 01014 
5.5  MeV 


3) ,EN(17) ,X(40) ,V(40) ,ANGLE(20) ,PI 


3), J  =  1,20) 
10.14431,-0. 
-11.43880,-0 
-12.73211,-0 
-14.38304,-0 
-16.22773,-0 
-17.55317,-0 
-18.67047,-0 
-18.90293,-0 
-18.69541,-0 
-18.68751,-0 

3), J  =  1,20) 
■10.06015,-0. 
-11.34385,-0 
-12.62155,-0 
-14.25599,-0 
-16.13103,-0 
-17.47239,-0 
-18.54096,-0 
-18.77677,-0 
-18.56277,-0 
-18.55547,-0 


/ 
97954,0 
.95306, 
.91973, 
.86783, 
.77307, 
.64425, 
.48302, 
.51104, 
.66378, 
.71336, 


00243, 
0.00245, 
0.00263, 
0.00329, 
0.00518, 
0.00752, 
0.00948, 
0.01016, 
0.01013, 
0.01010/ 


/ 

97856,0 

.95118, 
.91565, 
.86217, 
.76409, 
.63802, 
.49259, 
.52247, 
.67994, 
.73088, 


00253, 
0.00255, 
0.00273, 
0.00338, 
0.00523, 
0.00753, 
0.00947, 
0.01016, 
0.01013, 
0.01009/ 


3), J  =  1,20)  / 

9.96185,-0.97802,0 

-11.23784,-0.94879 

-12.48963,-0.91315 

-14.10896,-0.85694 

-16.01343,-0.75538 

-17.36967,-0.63302 

-18.39054,-0.50408 

-18.63231,-0.53571 

-18.41410,-0.69875 

-18.40473,-0.75195 

3), J  =  1,20)  / 

■9.87913,-0.96994,0 

-11.16730,-0.93379 

-12.41234,-0.89270 

-14.00646,-0.83481 

-15.92977,-0.73124 

-17.31061,-0.61073 

-18.35938,-0.48065 

-18.59678,-0.51740 

-18.31359,-0.70229 

-18.31906,-0.75424 


.00266, 

,0.00267, 

,0.00284, 

,0.00348, 

,0.00530, 

,0.00755, 

,0.00945, 

,0.01016, 

,0.01012, 

,0.01008/ 


.00284, 

,0.00285, 

,0.00303, 

,0.00365, 

,0.00542, 

,0.00763, 

,0.00952, 

,0.01026, 

,0.01017, 

,0.01013/ 
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DATA  ((PRIGAM(5,J,I),I  =  1 
C-8. 50918, -0.99050, 0.00311,- 
C-10. 37903, -0.95256, 0.00304 
C-l 1.04070,-0. 91026, 0.00312 
C-13. 01883, -0.85984, 0.00345 
C-14 . 76568 ,-0 . 78194 ,0 . 00450 
C-16. 51249, -0.66413, 0.00662 
C-17. 67830, -0.54866,0.00867 
C-18. 41481, -0.47923, 0.01010 
C-18. 18616,-0.64571,0.01026 
C-18. 06301, -0.77505, 0.01014 
4.5  MeV 

DATA  ((PRIGAM(6,J,I),I  =  1 
C-8. 35440, -0.99054, 0.00342,  - 
C-10. 20947, -0.95221, 0.00333 
C-ll. 45245, -0.90776, 0.00339 
C-12. 79416, -0.85529, 0.00371 
C-14. 51023, -0.77656, 0.00471 
C-16. 27127,-0.66035,0.00675 
C-17. 42659, -0.55583, 0.00871 
C-18. 11836,-0.50407,0.01011 
C-17. 92974, -0.66827, 0.01033 
C-17. 75403, -0.81884, 0.01016 
3.5  MeV 

DATA  ((PRIGAM(7,J,I),I  =  1 
C-8. 17183, -0.98828, 0.00389, ■ 
C-10. 03899, -0.94159, 0.00379 
C-ll. 28263, -0.88915, 0.00385 
C-12. 59264, -0.83133, 0.00414 
C-14 . 26622 ,-0 . 75042 , 0 . 00508 
C-16. 02860, -0.63648, 0.00700 
C-l 7 . 1 7563 , -0 . 54200 , 0 . 00887 
C-17. 85195, -0.50092, 0.01025 
C-17. 58929, -0.69379, 0.01046 
C-17. 32045, -0.87894, 0.01022 
2.5  MeV 

DATA  ((PRIGAM(8,J,I),I  =  1 
C-7. 92954, -0.98676, 0.00464, 
C-9. 79689, -0.93304, 0.00450, 
C-ll. 03843, -0.87135, 0.00454 
C-12. 31311, -0.80576, 0.00481 
C-13 . 91068 ,-0 . 72439 , 0 . 00565 
C-15. 62706, -0.61920, 0.00740 
C-16. 73064, -0.54145, 0.00913 
C-17. 37571, -0.51354, 0.01045 
C-17 . 13507 ,-0 . 70903 ,0 . 01078 
C-16 . 64530 ,-0 . 96460 ,0 . 01044 
1  5  MeV 

DATA  ((PRIGAM(9,J,I),I  =  1 
C-7. 58567, -0.98475, 0.00605, 


,3), J  =  1,20)  / 

-9 . 74493 , -0 . 96953 , 0 . 00305 , 
-11.02157,-0.93180,0.00306, 
-12.24270,-0.88884,0.00322, 
-13.80868,-0.82813,0.00382, 
-15.74444,-0.72261,0.00554, 
-17.14815,-0.60454,0.00769, 
-18.14754,-0.49304,0.00953, 
-18.37967,-0.53582,0.01029, 
-18.07745,-0.73343,0.01018, 
-18.06879,-0.79191,0.01012/ 

3), J  r.   1,20)  / 

9.57970,-0.96997,0.00335, 

-10.84418,-0.93049,0.00334, 

-12 . 04056 ,-0 . 88515 ,0 . 00349 , 

-13.56468,-0.82298,0.00406, 

-15.49461,-0.71747,0.00570, 

-16.90625,-0.60506,0.00777, 

-17.85917,-0.51282,0.00954, 

-18.11672,-0.55-524,0.01034, 

-17.78767,-0.76929,0.01022, 

-17.75121,-0.83932,0.01013/ 

3), J  =  1,20)  / 

9.40385,-0.96308,0.00381, 

-10.67691,-0.91572,0.00380, 

-11.86211,-0.86334,0.00394, 

-13.34028,-0.79744,0.00447, 

-15.24439,-0.69262,0.00601, 

-16.66515,-0.58478,0.00797, 

-17 . 59695 , -0 . 50551 , 0 . 00967 , 

-17.84378,-0.55745,0.01050, 

-17.38142,-0.81832,0.01031, 

-17 . 30779 ,-0 . 90388 ,0.01019/ 

3), J  =  1,20)  / 

•9.15743,-0.95848,0.00453, 

•10.43662,-0.90238,0.00450, 

-11.60502,-0.84194,0.00463, 

-13.02839,-0.76978,0.00511, 

-14.85561,-0.67139,0.00649, 

-16 . 23975 ,-0 . 57630 , 0 . 00829 , 

-17.12407,-0.51501,0.00987, 

-17.41302,-0.56089,0.01077, 

-16.80073,-0.87415,0.01058, 

-16 . 59248 ,-1 . 00247 , 0 . 01038/ 

,3), J  =  1,20)  / 

-8 . 81580 ,-0 . 94930 , 0 . 00590 , 
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C-9. 47027, -0.91570, 0.00586, -10. 12112, -0.87592, 0.00585, 
C-10. 72606, -0.83646, 0.00588, -11. 28246, -0.80051, 0.00594, 
C-l 1.95995, -0.75838, 0.00609, -12. 62827, -0.71931, 0.00634, 
C-13 . 43664 , -0 . 67503 , 0 . 00678 , -14 . 30309 , -0 . 62950 , 0 . 00748 , 
C-14. 99557, -0.59361, 0.00823, -15. 54272, -0.56631, 0.00899, 
C-15. 97945, -0.54493, 0.00971, -16. 32359, -0.53037, 0.01037, 
C-16. 58714, -0.52369, 0.01095, -16. 72687, -0.54101, 0.01139, 
C-16. 64413, -0.63341, 0.01169, -16. 23381, -0.81927, 0.01160, 
C-15. 88589, -0.96327, 0.01142, -15. 71597, -1.03426, 0.01132/ 

*  0.75  MeV 

DATA  ((PRIGAM(10,J,I),I  =  1,3), J  =  1,20)  / 
C-7. 14444, -0.99694, 0.00839, -8. 40557, -0.94672, 0.00815, 
C-9. 08480, -0.90034, 0.00808, -9. 7500 1,-0. 84825, 0.00805, 
C-10. 35795, -0.79770, 0.00806, -10. 90872, -0.75179, 0.00811, 
C-ll. 55655, -0.70042,0. 00822, -12. 17259, -0.65601, 0.00840, 
C-12. 87772, -0.61438, 0.00872, -13. 61158, -0.58111, 0.00921, 
C-14. 18642, -0.56337, 0.00974, -14. 63111, -0.55635, 0.01029, 
C-14. 97583, -0.55384, 0.01085, -15. 24170, -0.55479, 0.01136, 
C-15. 46310, -0.55269, 0.01186, -15. 64410, -0.55038, 0.01230, 
C-15. 85622, -0.54637, 0.01287, -16. 06595, -0.53718, 0.01350, 
C-16. 20228, -0.52879, 0.01393, -16. 26594, -0.52541, 0.01413/ 

*  0.325  MeV 

DATA  ((PRIGAM(11,J,I),I  =  1,3), J  =  1,20)  / 
C-6. 71592, -1.03172, 0.01153, -8. 07690, -0.94581, 0.01120, 
C-8. 78940, -0.88085, 0.01112, -9. 47427, -0.81217, 0.01108, 
C-10. 09263, -0.74702, 0.01108, -10. 63353, -0.69058, 0.01111, 
C-ll. 25494, -0.62881, 0.01120, -11. 82104, -0.57826, 0.01133, 
C-12. 43931, -0.53484, 0.01155, -13. 03722, -0.51037, 0.01186, 
C-13. 48116, -0.50724, 0.01218, -13. 80129, -0.51925, 0.01251, 
C-14. 03747, -0.53675, 0.01287, -14. 21578, -0.55331, 0.01323, 
C-14. 35568, -0.56531, 0.01359, -14. 46665, -0.57279, 0.01391, 
C-14. 60595, -0.57434, 0.01434, -14. 74812, -0.56676, 0.01483, 
C-14. 84392, -0.55716, 0.01517, -14. 89213, -0.55091, 0.01535/ 

*  0.1  MeV 

DATA  ((PRIGAM(12,J,I),I  =  1,3), J  =  1,20)  / 
C-6. 40158, -1.04235, 0.01675, -7. 81309, -0.92558, 0.01632, 
C-8. 50949, -0.85097, 0.01622, -9. 16255, -0.77588, 0.01617, 
C-9. 73855, -0.70714, 0.01617, -10. 23417, -0.64890, 0.01621, 
C-10. 78308, -0.58847, 0.01631, -11. 25874, -0.54288, 0.01643, 
C-ll. 74978, -0.50890, 0.01663, -12. 18135, -0.50036, 0.01687, 
C-12. 47044, -0.51439, 0.01710, -12. 66152, -0.54219, 0.01734, 
C-12. 78589, -0.57551, 0.01758, -12. 86513, -0.60814, 0.01784, 
C-12. 91613, -0.63496, 0.01810, -12. 95385, -0.65299, 0.01835, 
C-13. 00502, -0.66188, 0.01874, -13. 06711, -0.65253, 0.01923, 
C-13. 10413,-0.64177,0.01957,-13. 12132,-0.63549,0.01974/ 

*  energy  grid  for  tabulated  values 

DATA  (EN(I),I  =  1,17)  /. 10, .15, .2, .3, .4, .5, .6, .8,1. ,1.5,2. , 
C3. ,4. ,5. ,6. ,8. ,10./ 

*  ANGLE  GRID  FOR  TABULATED  VALUES 

DATA  (ANGLE(I) ,1=1,20)  /O. 5, 1.5, 2. 5, 4. 0,6. 0,8. 5, 12. 5, 17. 5, 25.0, 
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C35. 0,45. 0,55. 0,65. 0,75. 0,85. 0,95. 0,1 10. 0,130. 0,150. 0,170.0/ 
DATA  PI  /3. 14159265/ 

*  GUASSIAN  QUADRATURE  POINTS  FOR  GUASS  32  PT  QUAD 

*  ANGULAR  DIRECTIONS  X(I) 

DATA  (X(I) ,1=17,32)  /. 048307665687738, .144471961584796, 
C. 239287362252137, .331868602282128, .421351276130635, 
C. 506899908932229,. 587715757240762,. 663044266930215, 
C . 7321821 1874029 , . 7944837959679421 , . 84936761373257 , 
C . 896321 155766052 , . 93490607593774 , . 9647622555875059 , 
C. 98561 151 1545268,. 997263861849482/ 

*  GAUSSIAN  WEIGHTS. 

DATA  (V(I),I=17,32)  /9.654008851472801D-02, .095638720079275, 
C . 093844399080805 , . 091 173878695763 , 8 . 76520930044040 1D-02 , 
C . 08331 1924226947 , 7 . 81 938957870700 1D-02 , . 072345794108849 , 
C. 065822222776362, .058684093478536, .050998059262376, 
C . 042835898022227 , . 034273862913021 , . 025392065309262 , 
C .016274394730906 , . 00701861000947/ 

*  16  POINT  GAUSSIAN  QUADRATURE  DATA 

*  ANGULAR  DIRECTIONS  X(I) 

DATA  (X(I), 1=33, 40)/0. 095012509837, 0.281603550779, 0.458016777057. 
CO . 6 178762444026 , 0 . 755404408 , 0 . 865631 202388 , 0 . 94457023073 , 
CO. 9894009349916/ 

*  GAUSSIAN  WEIGHTS. 

DATA  (V(I), 1=33, 40)/0. 189450610455, 0.182603415045, 0.169156519395, 
CO. 14959598816,0. 124628971255,0.09515851168249,0.0622535239386, 
CO. 0271524594117/ 

END 
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CO-60  GAMMAS  USED  IN  BENCHMARK  EXPERIMTAL  SET  UP.   21  CM  CONCRETE  SHIELD. 

26   11   16  0.2099998E+02 

1.32910000  -0.21000000  0 . 0000000E+00  0 . 1250000E-02  0 . 254 6019E+00   1 
0.2500000E+02  0. 1000000E+04  0 . 5000000E+02 

1.24913500   1.08389700   0.94448430   0.83443190   0.72436540 

0.61940720   0.51417040   0.40399720   0.29370560   0.19372250 

0.09305161 

0.28694063E-01  0 . 12730098E+00  0 . 22590786E+00  0 . 27773422E+00  0 . 37065309E+00 

0.51540941E+00  0 . 67888516E+00  0 . 82364148E+00  0 . 91656035E+00  0 . 94387990E+00 

0.95959461E+00  0 . 98009801E+00  0 . 99581277E+00 

0.47022040E-05  0 . 15223846E-03  0 . 56299102E-03  0 . 95597235E-03  0 . 2052 1560E-02 

0.30115084E-02  0 . 52521937E-02  0 . 88381656E-02  0. 14977988E-01  0 . 17375331E-01 

0. 17090712E-01  0 . 48849441E-04  0 . 62234257E-03  0 . 88301837E-03  0 . 20534 184E-02 

0.28775800E-02  0 . 40205047E-02  0 . 67575760E-02  0 . 10567617E-01  0 . 16934209E-01 

0.19957155E-01  0 . 20271797E-01  0 . 20578112E-03  0 . 14748077E-02  0 . 21566362E-02 

0.26595094E-02  0 . 44032671E-02  0 . 52555390E-02  0 . 82526803E-02  0 . 12288887E-01 

0. 18553708E-01  0 . 22066455E-01  0 . 22938035E-01  0. 92780311E-03  0 . 22236416E-02 

0.28058987E-02  0 . 35442943E-02  0 . 49042068E-02  0 . 60939305E-02  0 . 91133118E-02 

0.13195600E-01  0 . 19330963E-01  0 . 23080699E-01  0 . 24239130E-01  0 . 44669174E-02 

0.48087537E-02  0 . 37399449E-02  0 . 54366104E-02  0. 64953342E-02  0 . 74979737E-02 

0.10756973E-01  0 . 14895651E-01  0 . 20586073E-01  0. 24771806E-01  0 . 26455171E-01 

0.18972654E-01  0 . 10200851E-01  0 . 73076747E-02  0. 83342344E-02  0 . 95182173E-02 

0.10106735E-01  0 . 13666026E-01  0 . 17567530E-01  0 . 22445098E-01  0 . 27324859E-01 

0.29720094E-01  0 . 43765388E-01  0 . 19039553E-01  0. 12353405E-01  0 . 12530528E-01 

0.13439607E-01  0 . 13594542E-01  0. 17267536E-01  0 . 20353712E-01  0 . 24325125E-01 

0.30082576E-01  0 . 33227455E-01  0. 66986382E-01  0 . 25518179E-01  0 . 17212339E-01 

0.17695744E-01  0 . 17719690E-01  0. 17160382E-01  0 . 20455401E-01  0 . 22344738E-01 

0.25482006E-01  0 . 32345820E-01  0 . 36221828E-01  0. 81144273E-01  0 . 28558638E-01 

0. 19103136E-01  0.21069471E-01  0. 21600537E-01  0. 20276256E-01  0 . 22555020E-01 

0.22941228E-01  0 . 26023451E-01  0 . 33709206E-01  0 . 38112491E-01  0 . 84768355E-01 

0.28907619E-01  0 . 21279372E-01  0. 20921603E-01  0. 22311699E-01  0. 21604430E-01 

0.23276970E-01  0 . 23180500E-01  0 . 26203763E-01  0. 34103289E-01  0 . 38666334E-01 

0.86422503E-01  0. 30600026E-01  0 . 20337187E-01  0. 21597095E-01  0 . 23381796E-01 

0.21648917E-01  0 . 23853570E-01  0. 23171041E-01  0 . 26310250E-01  0 . 34305081E-01 

0.38975243E-01  0 . 90028822E-01  0. 30385982E-01  0. 19603256E-01  0 . 22873476E-01 

0.24471726E-01  0 . 22968661E-01  0. 24191812E-01  0 . 22764482E-01  0. 26419207E-01 

0.34568090E-01  0 . 39366838E-01  0 . 91310024E-01  0. 29957492E-01  0 . 22471957E-01 

0.21923792E-01  0 . 24178460E-01  0 . 24028271E-01  0. 25184281E-01  0 . 22575960E-01 

0.26516844E-01  0 . 34756698E-01  0. 39668877E-01 

1   1.24500000   0.12130129E+00   1.00000000 
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ABSTRACT 

Radiation  that  is  scattered  in  the  air  back  to  a  ground  target  is  known  as 
skyshine  and  is  of  concern  in  any  radiation  facility.  For  unshielded 
gamma-photon  sources,  several  accurate  approximate  techniques  have  been 
developed  to  calculate  the  far-field  skyshine  doses.  However,  the  accuracy  of 
these  methods  when  applied  to  shielded  sources  remains  largely  unknown.  Thus, 
the  primary  purpose  of  this  work  is  to  develop  an  accurate  method  for  treating 
shielded  sources,  and  then  to  assess  the  accuracy  of  an  approximate  method  which 
uses  simple  buildup  and  exponential  attenuation  to  account  for  the  source  shield. 

The  resulting  composite  method  uses  an  accurate  one-dimensional  transport 
code  to  calculate  the  energy  and  angular  distribution  of  photons  emerging  from  a 
slab  shield  above  the  source.  Then  this  emergent  photon  distribution  is  used  as 
an  effective  point  source  by  an  approximate,  but  accurate,  method  based  on  beam 
response  functions  to  determine  the  far-field  skyshine  dose.  This  composite 
method  is  capable,  as  shown  from  comparisons  to  benchmark  experimental  data, 
of  accurately  calculating  the  shielded  skyshine  dose. 

The  composite  method  is  used  to  assess  the  accuracy  of  the  approximate 
MicroSkyshine  method  which  uses  buildup  and  exponential  attenuation  to 
account  for  the  source  shield.  Shield  thicknesses  from  0.001  to  6  mean-free-paths 
are  considered.  For  broadly  collimated  N-16  sources,  the  approximate  method 
results  are  within  a  factor  of  two  of  the  composite  results.  For  narrowly 
collimated  N-16  sources  the  approximate  method  is  generally  underpredictive  by 
more  than  a  factor  of  two.     At  source  energies  of  1.25  MeV  and  below,  the 


approximate  method  never  underestimates  by  more  than  a  factor  of  2,  and  tends 
to  overpredict,  often  by  a  factor  greater  than  two  (especially  for  broadly 
collimated  sources). 

The  composite  method  reveals  two  surprising  features.  First,  increasing  the 
thickness  of  a  thin  shield  (<  1  mean-free-path  length)  often  causes  the  skyshine 
dose  to  increase.  This  effect  is  especially  evident  for  high  energy  photons 
collimated  vertically  into  a  narrow  beam.  Second,  the  skyshine  dose  is  almost 
independent  of  source  photon  energies  for  source-to-detector  distances  less  than 
250  meters.  At  distances  greater  than  250  m,  the  skyshine  dose  increases  with 
increasing  photon  energy. 


